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Summary. — The k ■ p method is a semi-empirical approach which allows to 
extrapolate the band structure of materials from the knowledge of a restricted set 
of parameters evaluated in correspondence of a single point of the reciprocal space. 
In the first part of this review article we give a general description of this method, 
both in the case of homogeneous crystals (where we consider a formulation based 
on the standard perturbation theory, and Kane's approach) and in the case of non- 
periodic systems (where, following Luttinger and Kohn, we describe the single-band 
and multi-band envelope function method and its application to heterostructures) . 
The following part of our review is completely devoted to the application of the k ■ p 
method to graphene and graphene-related materials. Following Ando's approach, we 
show how the application of this method to graphene results in a description of its 
properties in terms of the Dirac equation. Then we find general expressions for the 
probability density and the probability current density in graphene and we compare 
this formulation with alternative existing representations. Finally, applying proper 
boundary conditions, we extend the treatment to carbon nanotubes and graphene 
nanoribbons, recovering their fundamental electronic properties. 
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1. Introduction 

To understand the physical properties of semiconductors it is necessary to know their 
electronic band structure, i.e. the behavior of energy as a function of the wave vector 
fc in the reciprocal lattice of the crystal. Several numerical methods can be successfully 
applied to find the band structure and the corresponding wave functions, such as the 
tight binding, the pseudopotcntial, the orthogonalized plane wave, the augmented plane 
wave, the Green's function and the cellular methods [1-3]. These methodologies can yield 
the desired results throughout the fc-space. 

Many phenomena, for example in the study of electrical transport (due to both elec- 
trons and holes) and of optical properties (such as absorption or gain due to electronic 
transitions caused by an incident optical wave), involve only the top of the valence band 
and the bottom of the conduction band. Indeed, low-energy electrons and holes arc sit- 
uated in these regions and also electronic transitions occur near the band edges of direct 
band gap semiconductors. Therefore a technique to obtain the band structure in such 
regions is of great interest. 

The k-p method [4-27] allows to extrapolate the band structure of materials from the 
knowledge of a restricted set of parameters (a limited number of energy gaps and of mo- 
mentum matrix elements between band lattice functions), evaluated in correspondence 
of a single point fcg of the reciprocal space, which are generally treated as fitting pa- 
rameters, that can be obtained from experiments or ab initio calculations. Even though, 
considering quite a large number of bands and thus of parameters, the fc • p method can 
be used to obtain the band structure all over the Brillouin zone of the material [28-32], 
its primary use is to explore with great detail the dispersion relations around the con- 
sidered point fcp. In particular, it allows to obtain the band structure of materials in 
the regions of the reciprocal space near the band extrema, expanding the eigenvalues 
and eigenvectors of the single-electron Hamiltonian as a function of fc around the wave 
vector fco corresponding to the band maximum or minimum. It has been shown to be 
very useful to study structures containing a large number of atoms, for which atomistic 
approaches would be computationally too expensive. 

This method, first introduced by J. Bardeen [33] and F. Seitz [34], was developed and 
adopted by W. Sliockley [35] and G. Dresselhaus, A. F. Kip and C. Kittel [36] in well- 
known papers on the energy band structures of semiconductors. It received a general 
formulation with E. O. Kane [10-12,37-39] and with J. M. Luttinger and W. Kohn [40,41]. 
It was later applied to strained materials (by G. E. Pikus and G. L. Bir [14]) and to 
heterostructures (for example by G. Bastard [42-44], M. AltarcUi [45-47] and M. G. 
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Burt [48-53]), proving to be a very useful and straightforward way to study the local 
properties of materials. 

In the last few years, a significant theoretical and experimental effort has been de- 
voted to the study of graphene and graphene-related materials, which appear as promis- 
ing candidates for many technological applications and are characterized by very un- 
usual and interesting physical properties. In particular, the application of the k ■ p 
method to the study of the electronic properties of graphene, sistematically pursued by 
T. Ando [54-56] and other authors, results in a description of the graphene properties 
in terms of the Dirac equation, the same relation that describes the relativistic behavior 
of elementary spin-(l/2) particles. This is at the basis of the experiments aiming to 
observe in graphene, at non-relativistic speeds, the analogous of purely relativistic quan- 
tum phenomena [57-61]. The application of proper boundary conditions to the relations 
found for a sheet of graphene allows to obtain the electronic properties of carbon nan- 
otubes and graphene nanoribbons, materials which (contrary to unconfined graphene) 
can exibit (depending on their geometrical details) a non-zero energy gap. 

The first part of this review is a short introduction to the k ■ p method in some of its 
most common formulations. 

In particular, sect. 2 describes the application of the fc • p method to homogeneous 
crystals, where, due to the periodicity of the material, the electron wave functions are 
Bloch functions and thus the unperturbed Bloch lattice functions arc adopted as a basis 
for the method. We first describe (following W. T. Wenckebach [5]) how the k p approach 
can be derived by just applying the standard perturbation theory to the Schrodingcr- 
Bloch equation and how this formulation can be adopted to study the dispersion relations 
of semiconductors with the diamond or zincblcndc structure. Then we briefly summa- 
rize the alternative formulation by Kane, consisting in the exact diagonalization of the 
Schrodinger-Bloch Hamiltonian for a subset of bands, and in the inclusion of the effect 
of the other energy bands with the Lowdin perturbation theory. 

In sect. 3, instead, we describe how the k ■ p method can be applied to the case of 
non-periodic systems, where the phase factor (involving the wave vector measured from 
the considered extremum point) of the Bloch lattice functions has to be replaced by 
properly defined envelope functions. Following J. M. Luttinger and W. Kohn, we derive 
the single-band and multi-band envelope function equations, and then we briefly outline 
the main approaches followed in the application of the envelope function theory to the 
study of heterostructures. 

The second part of the review is devoted to the application of the k ■ p method, 
and in particular of the envelope function approach, to graphene, carbon nanotubes and 
graphene nanoribbons. 

In sect. 4, following T. Ando, we perform a first-order expansion of a simple tight- 
binding description of graphene, obtaining the Dirac equation for the envelope functions 
(corresponding to the two degeneration points of graphene) in the presence of a generic 
external potential, and we analytically solve this equation for the case of null potential. 
Starting from this formulation, we also derive general expressions for the probability den- 
sity and for the probability current density in graphene, and we compare them with those 
used, adopting a valley-isotropic representation, by C. W. J. Beenakkcr ct al. [62,61]. 

In sect. 5 we extend the previous treatment to the study of carbon nanotubes, en- 
forcing a periodic boundary condition along the chiral vector, that univocally charac- 
terizes these tubules. In particular, we show how this periodic condition on the overall 
wave function translates in terms of the envelope functions, and we analytically solve the 
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Dirac problem in the absence of an external potential, obtaining the conditions for which 
nanotubes have a semiconducting or a metallic behavior. 

Finally, in sect. 6 we discuss the case of graphene nanoribbons. Adapting the approach 
adopted by L. Brey and H. A. Fertig [63,64] to the mathematical formulation of graphene 
proposed by T. Ando, we study both zigzag and armchair nanoribbons, obtaining an 
analytical solution in the absence of an external potential, and recovering the fundamental 
properties of these structures. 

2. The k p method in homogeneous crystals: derivation based on the stan- 
dard perturbation theory and Kane's model 

We begin our overview of the k ■ p method describing its formulation in the case of 
homogeneous crystals. 

In a pure crystal an electron is subject to a periodic potential energy 

(1) ULir)^ULir + R), 

with R any linear combination of the lattice vectors, and thus also the Hamiltonian 
is invariant under translation by the lattice vectors. Therefore, if ip^ir) is the wave 
function of an electron moving in the crystal, also tp^{r + R) will be a solution of 
the Schrodinger equation and therefore will coincide with ip^]^{r), apart from a constant 
with unit modulus (otherwise the wave function could grow to infinity, if we repeated 
the translation R indefinitely). Thus the general form of the electron wave functions 
will be 

(2) i^Ur)=e^''-ul{r), 

where "ip^ir) is usually called "Bloch function", while w^(r) is called "Bloch lattice func- 
tion" and is periodic with the lattice periodicity 

(3) ul{r + R)^uUr) 

(Bloch's theorem) [65]. 

Starting from the Schrodinger equation (in the absence of a magnetic field) for ip'j^{r) 

(4) H^°UUr) = E^i'Ur), 
with (in the absence of a magnetic field) 

(5) i7(0) = _^V2 + C/i(r) 

(where rrie is the electron mass and h is the reduced Planck constant) and substituting 
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■0^'(r) with the generic expression of the Bloch function, we obtain 

(6) f-^v^ + uUr))e^''-uUr) = 

V 2 me / 



■(e^'=-V.(V<(r)+^fc<(r)) 



2 rUe 

+ (Ve^'=-) • (V<(r) + ikul{r))) + UL{r)e^^-''ul{r) - 



(e^^-%V^ul{r)+zk- Vulir)) 



2 me 

+ (zfee^'=-) . (V«J!(r) +zfc<(r))) +C/i(r)e''^-<(r) 



e"=-(v2zij!(r) + zk ■ Vulir) + ik ■ VuHr) - k^uUr)) 



2 m 

+ ULir)e^'^-^ulir) = 



V' + ULir) k-V + ^ ulir) 



2 rrie J me 2 



and thus 

(7) Hulir) = (ijW + i/(i))uj!(r) = E^^uHr), 
with 

(8) i/(i) = _l^fe.V 



me 2 Trie 

(where k = |fc|). What we have just obtained is clearly an equation for the Bloch lattice 
functions (the Schrodingcr-Bloch equation), which needs to be solved only for a single 
primitive cell with the boundary condition that the function uHr) must be periodic 
with the lattice periodicity. For each value of k this equation has a periodic solution 
only for selected values EJ^ of the energy E. Noting that 7J(^)(r) reduces to zero when k 
approaches and thus that this part of the Hamiltonian can be treated as a perturbation 
around A; = 0, we can locally solve this equation using the time- independent perturbation 
theory, assuming to know the eigcnfunctions and eigenvalues of H^^\r), i.e. the Bloch 
lattice functions and the energy band values for fc = 0. 

For most of the semiconductors the maximum of the valence band is in the F-point 
(the center of the first Brillouin zone represented with the Wigner-Seitz method) and 
therefore corresponds to fc = 0; the minimum of the conduction band instead is for 
fc = only for direct-gap semiconductors. When the extremum point of the energy band 
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(and thus the interesting region) is for a generic ko, we can easily extend this argument 
observing that, if we define the value of H in feo as 



(9) 



Hio) 



-fco • V 



2mp 



we have that the value of H in k is 



(10) H = H^"^ + H^'^ = Hk„ 



Hko 



(fc-fco)-V 



■(fc-fco)- V- 



h'{P-ko') 

2 TTlp 



2 mp 



— (fc — ko)-hko (fc — ko)-hko 



me 
h ^ 
me 

h ^ 

rrie 
h ^ 



ko) 


{hko 


-ihV) + 


ko) 


{hko 


-inv) + 


ko) 


(hko 


-^inv) + 


ko) 


{hko 


-inv) + 



2 TTlt 
2 7Tl( 
2 TOj 
2 TO, 



ko^ - 2fc • feo + 2fco^) 



■(fc2-2fe-fco + A:o^) 
■(fc - fco) ■ {k-ko 



|fc-fco|' 



and for fc near fco the term between square brackets can be treated as a perturbation of 
Hko [^^]- sake of simplicity, in the following we will consider feo = 0. 

An important point to notice is that, for any selected fc, the functions u^{r) form an 
orthogonal and complete set (in the restricted sense that any function with the lattice 
periodicity can be expanded in terms of the Bloch lattice functions corresponding to the 
selected fc). 

To describe the main results of time- independent perturbation theory [66,67], we have 
to distinguish the case in which the unperturbed energy levels are non-degenerate from 
the case in which such a degeneration exists (in the following we will use the notation 
of W. T. Wenckebach [66]). Let us begin from the first case. The problem we have to 
solve is 



(11) 



where H^--^'^ is the unperturbed Hamiltonian and _ff the perturbation. If we expand the 
eigenvalues £"„ and the eigenfunctions \n): 



(12) 



|n) = + + \n)^'^'> + 



we insert these expressions into the eigenvalue equation, and we enforce the identity 
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between terms of the same order, we find 

(13) = 

i?(0)|7i)W+i7(i)|n)(°) = i?(0)|n)«+i?W|n)(0), 
i7(")|7i)(2)+ij(i)|n)(i) = i?(0)|n)(2)+ij(i)|n)(i)+i?(2)|„)(0)^ 



The first equation corresponds to the unperturbed eigenvalue equation, the solutions of 
which, En'' = Eq and = |7i0), are assumed to be known. From the other equations, 

instead, we can obtain the corrections to these values produced by the perturbation H^^\ 
In particular, if we stop to the first-order corrections for the eigenfunctions and to the 
second-order corrections for the eigenvalues we find 

(14) \n) ~ |7i0) + |7i)(i) 

(choosing {nO\n)^^'' = 0) and 

(15) En ^ EJi + El'^ + Ei'^ 



(TOO|ff(iVo) 

rpn Tprn 



^ E^ + {nO\H'^^^\nO) 



/(nO|7?(i)|TOO)(mO|i?(i)|nO) 

rn^n 



E"-E" 



When we examine degenerate unperturbed states, the expressions we have just found 
diverge and thus we have to modify our treatment. In particular, if the degenerate energy 
level Eq corresponds to a multiplet of degenerate states \naO) (with a = 1,2,...,(7„, 
where gn is the degeneracy) and we have to solve the perturbed problem 

(16) H\^) = [H^°^ +H^^^M =E\i;), 
we can express the new generic eigenfunction 1-0) as 

(17) IV-) =^|na)(na|^), 

a=l 

where the \na)'s are states which are related to the unperturbed eigenvectors |naO)'s by 
the perturbation matrix elements between different multiplets (as we will see in eq. (23)). 
If we define 

(18) Hl\ = {7ia\H\nb) {7m\[H^^'^ + H^^'Unfe), 
we can express our perturbed equation in the following way: 

(19) Y.^ab{nbm=E{na\^). 

b=l 
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Noting that the definition of the i?"f,s can be equivalently expressed in this way 

(20) [if(")+i/(i)]|n6)=|i|na)if,V 

a=l 

inserting into this equation the expansions 

(21) h:, = (ij:j(«) + (i?:j(i) + (ij,"j(2) + . . . , 

\na) = |na)(o) + |na)(i) + |na)(2) + ..., 
and enforcing tlie identity of the terms of the same order, we find 

(22) 77(°)|n6)(")=.£|na)(°)(77:,)(°), 

a— 1 a— 1 

+ = 1^ |na) (//^",)(o) 

+ £ |na)W(i7„"jW + £ |na)(°)(i7;'j(^), 



The first equation corresponds, noting that (i^^fc)'''^'' ^ E^Sab, to the unperturbed eigen- 
value equation, the solutions of which, Eq and \na)^^'' = \naO), are assumed to be 
known. From the other equations, instead, we can obtain the corrections to these values 
produced by the perturbation. In particular, if we stop to the first-order corrections for 
the eigenstates and to the second-order corrections for the eigenvalues, we find 

(23) \nb >^ \nbO) + \nb)(^^ = \nbO) + J] £ (|mcO) ^''^ ) 

m^n c=l ^ ^0 - ^0 / 

(choosing (ncO|n6)^^) = 0) and 

(24) H:^, c (iJ,"j(«) + (i75,)(i) + (i?,"j(2) ^ j^ng^^ ^ („cO|ij(i) \nbQ) 



EE 

ni^n a—1 



Once the iJ"^ have been found, we can obtain the energy levels E solving the equation 



(25) 



Y,H:,{nb\i')=E{na\i'), 

6=1 
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or, equivalently, find-ing the eig6nva,lues of the nicitrix U^^ (matrix x (j^i with elements 
Hab) by solving 

(26) det (iJ" -£;/)= 

(with / the gn x gn unit matrix). We notice that, computing also the eigenvectors {na\ip) 
of such a matrix and combining such results with the \nb) that have been computed 
before up to the first order, it is also possible to know the eigenfunctions of the 
perturbed problem. 

In the case of the k ■ p Hamiltonian that we have found before [5] , we can use the 
UQ(r) (ufc(r) for A; = 0) as |nO) and we have that 



(27) {mO\H'-^^\nO) = {mO\ 



■ t;2 

fc- V 



2 7,2 



h k 

hiO) + (mOl- ItiO) = 

2 ITLp 



— ■ {mQ\{~ihV)\nO) + (mO|^-^|nO). 

TOe 2 me 

The second term clearly gives only diagonal matrix elements, because it is equal to 
{h'^k'^/(2me))6nm- The first term, instead, gives only non-diagonal matrix elements 
because it is known [68] that 

(28) {nko\{~ihV)\nko) + hko = ra^Vn = -^VkE^ 

(where Vn is the expectation value of the velocity of the Bloch waves, and in our consid- 
erations we are assuming fep = 0) and Vfc£'^* = in the band extrema. 

Then, if the unperturbed energy bands are non-degenerate, we can write that 

p„ p„ ^ ^ (nO|fc■(-^?iV)|77^0)(mO|fc■(-^?lV)|nO) 



rn^n 

2 ^-^ m* 



^0 



-fin 



where fi^v = x, y, z, while m* is the effective-mass tensor defined by 



-1-1 r\ pnm pmn 

(30) ^ = —'5.. + A E 



and the momentum matrix elements at the band extremum are 
(31) P-^"^{nQ\{^^nW,)\mO). 
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If the unperturbed energy bands are degenerate, instead, we have 

(32) {Hi:)eb = E^Scb + t^Seb + —{ncQ\k ■ {^tnV)\nbO) 

^ ^ {ncO\k ■ {-ihV)\maO){maO\k ■ {-ihV)\nbO) 



where fijiy = x, y, z, while mj^^ is the effective- mass tensor defined by 

^ rn^n a—1 ^ ^ 

and the momentum matrix elements at the band extremum are 



(34) (P^)™^ - {ncO\i-ihV^)\mbQ). 

In most of the cases all the {Pf_i)cb — 0, and the linear term in kp, disappears. The energy 
levels will be found solving 



(35) det {HI - EI) = 0. 

Thus, in principle to perform a calculation of the energy bands we would have to know 
the |nO)'s (the Bloch lattice functions at fe = 0). Since the Hamiltonian H'^^^ and 
its eigenfunctions |nO) have the periodicity of the lattice, the problem can be solved 
inside a single primitive cell, enforcing periodic boundary conditions at the surface of the 
cell. Most semiconductors of interest have the diamond or zincblende crystal structure; 
for these materials we can choose as lattice primitive cell a Wigner-Seitz cell centered 
around an atomic site (the one with the strongest potential in the case of the zincblende 
structure, characterized by atoms that are not all identical) and with, at four vertices 
of the cell, four other atoms forming a tetrahedron with the center coincident with the 
primitive cell center (fig. 1). We can use a central force model (the same results can be 
obtained using group theory), considering the potential inside the primitive cell as due 
only to the attraction of the nucleus of the central atom, shielded by its electrons [5]. We 
find that the Bloch lattice functions at fc = exhibit symmetry properties analogous to 
those of atomic orbitals: we have completely symmetric s-type states pys(r), and p-type 
states antisymmetric with respect to a coordinate and symmetric with respect to the 
others, i.e. of the form p,jx{r)x, Pijy{r)y, and p^z{r)z (where r = i/aJ^ + y^ + z^) . Then, 
treating the electrostatic potential of the cores at the vertices of the primitive cell as a 
perturbation, we see that, to first order, this potential does not change the eigenfunctions 
but shifts the energy levels and in particular breaks the degeneracy between each s-type 
state and the corresponding three p-type states (which remain mutually degenerate) . As 
a result, we find that at A; = the top of the valence band can be described with three 
degenerate states: \vxO) = pv(r)x, \vyO) ~ Pvi'r)y and \vzO) = p^{r)z, while in most 
cases the bottom of the conductance band is described by a non-degenerate symmetric 




state |cO) = Pc{r) (with the important exception of silicon, where at fc = also the 
bottom of the conduction band is characterized by three states \cxO), \cyO) and \czO)). 
Therefore, if we treat the conduction band as a non-degenerate band, we obtain 



1% X ^ k^kjj 



(36) K = E'o + tT. 



where /i, = x, y, z and 

1 _ 1 , 2 ^ (cO|(-^?iV^)|mO)(mO|(-^?^V,)|cO) 



me 



The largest contribution to the sum comes from the bands m for which \E^ — 
is smallest, i.e. from the three valence bands. If we compute the momentum matrix 
elements between the valence bands and the conduction band, we find that, due to the 
symmetry properties of the Bloch lattice functions. 



(38) (w/iO|(-inV^)|cO) = -(cO|(-i?iV^)|u/^0) = -ihPS, 



with ^,i> = x,y, z and P = (u/iO|Vp|cO) a non-zero quantity, which multiplied by h has 
the dimensions of a momentum. Consequently, the effective mass in the conduction band 
that we find is isotropic and equal to 



with E°=E^-E^. 

As to the valence band, we must use the degenerate perturbation theory and, with a 
motivation analogous to that used in the study of the conduction band, we can consider 
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only the interaction between the three degenerate valence bands and the conduction 
band, which is the nearest energy band. Thus, using the previous results, we have that 

(40) (^^^W = i?o'5«, + f 



with 



1 1 , , 2 ^4p.{vaQ\{-ihV,)\maQ){maQ\{-ihV^)\vpQ) 
(41) ^ = — -^"/J-^M^ + ZTTZ^ 



^0 ^0 



1 ^ 2 (waO|(-i?iV^)|cO)(cO|(-i?iV^)|?;/30) 

-OapO^i, H J 



m, 
1 

m, 



^0 -^0 



— da/Sd^^ - ^ 2^0 ^"A''5/3>^ 



and thus the valence energy bands near the extremum can be obtained finding the eigen- 
values of the matrix 



(42) 



21.2^ 



HI 



2 me 



kx ky kx 



Till now we have not considered the effect of the so-called spin-orbit interaction, which 
often has a non-negligible influence on the energy bands. The physical phenomenon is 
the following [69,70]. An electron has an intrinsic magnetic moment 

h h eh <T 

(43) M = -7e- (T = ~g^j^-(T = ~ge- -a- = -gefJ-B 77' 

Z 2 2 rrie 2 2 

where e is the modulus of the electron charge, cr is a vector with three components 
consisting of the Pauli spin matrices: 

(44) = ( ? J ) ' ^y-{] V ) ' = ( -1 ) ' 

7e is the intrinsic gyromagnctic ratio of the electron, 7^ is its orbital gyromagnetic ratio, 
<7e = 7e/7L is its intrinsic g-factor and = eh/(2me) is the Bohr magneton. When 
an electron moves in a system (such as the atom) where the charge distribution (for 
example the nucleus charge) produces an electric field E, for the theory of relativity this 
electric field will appear as a magnetic field in the frame of reference of the electron. In 
particular if the motion of the electron were uniform the equivalent magnetic field would 
be equal to B = —{v x E) /c^ . The fact that the electron (and its frame of reference) is 
rotating halves such an equivalent magnetic field [69,70]. Thus the Hamiltonian of the 
electron will have an additional part 

/.-N / E X V\ eh , h ,,„-rr N 

(45) ffso = MB • — ^ = -(Exv)^ ■ {{VUl) x v) 

\ 2c J 4 nieC 4 irieC 
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(with U L the potential energy) , which in the absence of an external magnetic field can 
be written also as 

(46) Hso = -r\^T ■ {{VUl) x p). 

However, if we insert this additional term into the original Schrodinger equation for the 
wave function ip^^{r) = e'''''^u^{r), we obtain 

(47) Hsomr) = ■ H^Ul) x ^hV)) e^^W^ir) = 

h 



4 m^c^ 



ik-r 



,,_ h 



a ■ ((VC/i) X fe) + ——^rr ■ ((VC/^) x {-ihV)) ul{r). 



If we repeat the procedure used to move from the Schrodinger equation for the wave 
functions i!^{r) to the Schrodinger-Bloch equation for the Bloch lattice functions uj!(r), 
we obtain that in the Hamiltonian of this last equation there will be two additional terms: 

(48) -^<T . {{VUl) xk) + ■ UVUl) x {-z h V)) = 
-—^a-{{VUL)xk)+Hso- 

The first term near fc = is small compared with the other term; thus only the second 
term is usually considered. The second term in the case of a potential energy with 
(locally) spherical symmetry (and thus of a radial electric field) becomes 

eh , „ , eh Er , 

(49) Hso ^ -^^(T ■ {E X p) ^ —^(T ■ ^{r X p) = 

eh^Er \ , .A , _, 

cr • (r X V) = —i—(T ■ (r x V). 



4 m'ic^r I 2 



In order to calculate the influence that the spin-orbit term has on the valence bands, we 
need to calculate the matrix elements on the basis states |wxO), |wyO), \vzQ) and |cO). 
Due to the symmetry proprieties of such states, we see that the only non-zero elements 
are the non-diagonal elements between valence band states 

(50) {vyQ\Hso\vxQ) = ~{vxQ\Hso\vyQ) ^iXa,, 
{vzO\Hso\vyO) = -{vyO\Hso\vzO) = iXa^, 
{vxO\Hso\vzO) = -{vzO\Hso\vxO) =iXay, 

with A a non-zero quantity given by (if Vc is the volume of the lattice unit cell) 

(51) x=^yj^x^pl{r)dr. 
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Therefore, considering also the spin-orbit coupHng, the matrix Hj^ becomes 



(52) 



TTV 



+iX 





]'- 













ay 







— Ox 










kx ky kx k^ 



fr2 



ky 



where ax, ay and tr^ are the Pauh spin matrices (44), which do not commute with one 
another. If we consider the special case fe |j i we can quite easily find the eigenvalues 
of this matrix, arriving to a third-order equation in the energy, the solutions of which 
represent the dispersion relations of the three valence bands, each one degenerate with 
respect to the spin. In particular, if we make the approximation {h*P^k'^ /{irilEg)) ^ A, 
we find the solutions [5] 



(53) 



Ehh = Eq 



X + — — k^, 

2 nip 



° 2 me V 3m,E0 

2h^ P2 



E 



Xh 



2A 



2 me 



3 rueEO 



Thus, considering the effect of the spin-orbit interaction, we have obtained (fig. 2) two 
valence bands (the heavy-hole band and the light-hole band) degenerate at fc = 0, where 
they have an energy Eg — E^ — (£'^ -I- A) = Eg — X lower than the conduction band, and 
one valence band (the spin-orbit band) which for k = has an energy A = 3 A lower than 
the other two valence bands. We notice that, while the light-hole band and the spin- 
orbit band have a negative effective mass of the same order of magnitude as the effective 
mass of the electrons in the conduction band, the heavy-hole band is characterized by a 
much larger effective mass (the fact that the obtained effective mass is positive instead 
disappears with a more refined treatment: obviously the effective mass of the electrons 
in the valence bands has to be negative, which corresponds to a positive effective mass 
for the holes). 

This simplified model is amenable to several refinements. 

As to the conduction band, we can include in the calculation the spin-orbit splitting 
of the valence band and the effect of the higher conduction bands. In particular, with 
the first change we obtain a better expression for the effective mass in the conduction 
band: 



(54) 



1 

m* 



1 



h^P^ 



3 Eg 'S{Eg+A) 



where Eg = E'^ - X. 
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Fig. 2. - Band structure near fe = obtained with the simplified model described in the text 
(the heavy- hole band has a wrong curvature) (adapted fiom [5]). 



Also in the treatment of the valence bands we can consider the effect of the higher 
conduction bands; one of the effects is that the resulting valence bands lose their isotropy 
and exhibit a complex orientation dependence in the reciprocal space ("band warping"). 

It is important to notice that the expressions found for the band structure depend on 
a small number of parameters, for example Eg, A and m* (from which we can calculate 
the parameter P using the expression found for the effective mass of the conduction 
band). From a practical point of view, these quantities are commonly obtained from a 
priori band structure calculations or, better, experimentally: in particular the bandgap 
values Eg and A are accurately known from optical experiments, while m* is known from 
cyclotron resonance experiments. 

The approach based on the "traditional" perturbation theory, that we have reported 
in this first part following the description of T. Wenckebach [5], differs from the method 
proposed by E. O. Kane [10-12,37-39]. 

Starting from the consideration that the Bloch lattice functions can be expanded 
in terms of the complete, infinite set of the unperturbed Bloch lattice functions, Kane 
computes this expansion in an approximate way, considering only a finite set of bands. 
In particular he considers only the three valence bands and the conduction band (not 
including the effects of the other bands) and diagonalizes exactly the Hamiltonian of the 
Schrodinger-Bloch equation in the presence of spin-orbit interaction [38] , written taking 
as a basis the following set, made up of a linear combination with constant coefficients of 
the Mo(r) considered in the absence of spin-orbit {i.e. of the functions |cO), |i;a;0), \vyQ) 
and \vzO) taken with spin-up and spin-down): 



(55) z|cO)|;), l/V2i\vx0)-i\vy0m), \vzO)\i), -1/V2{\vx0) + i\vyO))\t), 
i|cO)|t), ~l/V2i\vxO)+i\vyO))\i), \vzO)\t), 1/V2i\vx0) ~ i\vyO))\l) 



(where [f) and W) are, respectively, the spin- up and spin-down unit spinors). 

From this diagonalization he finds, for small values of fc^, the following expressions 
for the considered bands (choosing the zero of energy at the top of the light-hole and 
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heavy- hole bands and defining tlie various quantities as before): 



(56) 



Ehh 



Eih 



E,_ 



'C 




+ 



3me(£;g + A) 



) 




Exh 



-A + 1 

2 TOe V 



3TOe(-B<,+A) 



2h^ P2 



) 



These expressions are very similar to the expressions obtained with the previously de- 
scribed simplified model, but clearly show the dual effect that each reciprocal interaction 
has on the related couple of bands. As before, these results give an incorrect effective 
mass for the heavy-hole band. 

From the diagonalization Kane also finds the Bloch lattice functions M^('r) that diag- 
onalize the Hamiltonian of the Schrodinger-Bloch equation in the presence of spin-orbit 
interaction {i.e. the eigenfunctions of this Hamiltonian) as linear combinations of the 
Mg (r) considered in the absence of spin-orbit; in particular for vanishing k they are (in 
the simplest case in which k \\ z): 



-l/V2{\vxO)+i\vyO))\t), l/V2{\vxO)-i\vyO))\i), 
l/VQi\vxO) - i |«2/0))|t) + v/273|wz0)|i), 
-1/V6(|^;x0) +1 |z;yO))|;) + ^fzJlWm). 
l/^/3(|ra0) -z|i;yO))|t) -l/\/3|i;z0)|i), 
1/V3(|wa;0) +z|?;2/0))|i) + l/y3|i;z0)|t). 



In order to take into account the effect of higher and lower bands on the considered 
ones, Kane uses the Lowdin perturbation theory [71,72]. Following this method, one can 
divide all the bands into two sets A and B: A is the set we want to treat exactly and 
_B contains all the other bands. At the lowest order of perturbation theory the coupling 
between the set A and the set B can be removed introducing the perturbed functions 



(57) 



i|cO)|i), 



*|cO)|t), 



(58) 




where z is in ^ and n is in _B. The renormalized interactions connecting m' 
given by 



and u[: are 



(59) 
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(with i, j in A). In this way we can reduce the Hamihonian matrix, which in principle 
connects aU the possible bands, to a Hamiltonian matrix relating only the bands of inter- 
est, but in which, however, the interactions with the non-considered bands are included. 
The method is accurate as long as <?C {Ha — Hnn\, with i in A and n in B, and 

thus the set A has to be selected in order to satisfy this relation (for example, also states 
degenerate with those in which we are interested have to be considered inside the set A). 
Note that the Lowdin perturbation theory reduces to the ordinary perturbation theory 
when only a single band is considered in the set A. 

Kane applies this perturbation method, starting from the Bloch lattice functions (57) 
of the set A of considered conduction and valence bands and from the unperturbed 
Bloch lattice functions of the set B of the higher and lower bands, obtaining a better 
approximation of the actual dispersion relations of the considered bands. 

An exact diagonalization of the Hamiltonian has also been performed (originally by 
M. Cardona and F. H. PoUak [28], more recently by other authors [29-32]), extending the 
number of considered bands (and thus the number of involved parameters) to reproduce 
the band structure all over the Brillouin zone to a reasonable degree of accuracy (for 
example, in their original paper M. Cardona and F. H. Pollak consider 15 bands, with 
10 parameters, to reproduce the energy band structure of germanium and silicon). 

3. The k ■ p method in non-periodic systems: envelope function theory and 
application to heterostructures 

Till now, we have considered the applications of the k ■ p method to periodic, ho- 
mogeneous crystals, using a description of the electron wave function in terms of Bloch 
functions. However, in the presence of a generic external potential the periodicity of 
the potential inside the crystal breaks down and thus the electron wave functions are 
far from periodic. Since the Bloch functions \nk) = e''° ''u^(r)/A/(27r)^, considered as 
functions of r and fc, are a complete set of orthonormal functions, also in this case the 
generic wave function could be expanded on the basis of Bloch functions in this way 



(where the sum over the number of bands together with the integral over the Brillouin 
zone corresponds to an integral over all the reciprocal space). However, in general a 
large number of Bloch functions, evaluated over a large range of wave vectors, would be 
necessary in this expansion. Therefore in this case it is convenient to replace the Bloch 
phase factor, involving the wave vector measured from the reference extremum point, 
with an envelope function, and thus to use a different formulation of the k ■ p method, 
based on the concept of envelope functions (^). 

In order to introduce this concept, we can make a very approximate calculation [77,78] 
in the hypothesis that the external potential energy U{r) ("external" here meaning "not 
due to the periodic structure of the lattice") is slowly varying on an atomic scale and 
the n-th energy band that we are considering is non-degenerate (thus with unique in- 
dependent Bloch lattice function u^{r)). In this case, the Schrodinger equation (in the 



(^) Notice that there is also an alternative approach to the envelope function theory using the 
definition of Wannier [73] and Slater [74], based on Wannier orbitals. See also [52,75,76]. 



(60) 




n 
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absence of a magnetic field) for the electron wave function ^'i^) 

(61) (-TT—"^^ + UL(r)] i'ir) + U(rU(r) = H^°^{r) + U(rU(r) = Ei:(r) 
\ 2me J 

(where UL{r) is the periodic lattice potential energy and ij'^-* is the Hamiltonian in the 
absence of the external potential energy U{r)) is equivalent to the equation 

(62) En{~iV)F{r) + U{r)F{r) = EF{r), 

where En{~i'^) represents the operator obtained replacing, in the dispersion relation 
En{k) describing the n-th energy band in the absence of the external potential, each 
component of k with the corresponding component of —i V, and F{r) is the envelope 
function, a slowly varying function that, when we consider only the rt-th band, multiplied 
by the fast varying Bloch lattice function Uq (t) (considered in fc = 0) gives the electron 
wave function. 

Indeed, if we expand i'ir) in the orthogonal basis set \i^k) = e^'''^u'(,(r)/VV (with V 
the crystal volume) 

(63) V'(r-)=^a,(fe)|i.fc), 

we can re-write the Schrodinger equation (61) in matrix form using the basis ji^fc) 

(64) J2 ((^^^1^*°^ + C/(r)|i^'fc>^'(fe')) = Ea,,{k) ^ 
E,{k)a,{k) + {{vk\U{r)\v'k')a„,{k')) = Ea^k), 

u' .k' 

where we have used the fact that (being \vk) an eigenfunction of H^^^ with eigenvalue 
EAk)) 

(65) (i/fc|ij(0Vfe') = E,,{k'){vk\v'k') = E,{k)5,,.5kk'- 
In particular, for v = n we have that 

(66) En{k)anik) + {{nk\U{r)\,^'k')a,,ik')) = Eanik). 

u' .k' 

If instead wc expand the envelope function equation in the orthogonal set of plane waves 
|A;) = /VV 



(67) 



i^(r) = ^a(fe)|fc), 

k 
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we can re- write the envelope function equation (62) in matrix form using the basis |fe) 

(68) ^ {{k\En{-i V) + U{r)\k')a{k')) = Ea{k) 

k' 

E„ik)aik) + J2 {{k\U{r)\k')aik')) = Ea{k), 

k' 

using the fact that 

(69) i-iv.r\k') = (-iv,)f(e*^'-7%/y) = (kirie^'^'-yVv) = {KY\k'), 

witli V — x^y, z and thus 

(70) En{-iV)\k') = E.n{k')\k') 

(being En{—i^) an operator made up of operators of the type {—iV„y) and then 
exploiting the orthogonality relation {k\k') = S^k'- The two equations (66) and (68), 
obtained from the Schrodinger equation and from the envelope function equation are 
exactly equal if 

(71) {nk\U{r)\,y'k') ^ 5n,'{k\U{r)\k'), 

i.e. if the matrix elements of the external potential U(r) between states from different 
bands are negligible. This is what happens if U is slowly varying on an atomic scale. 
Indeed, in this case we have that 

N 



(72) (nfc|C/(r)|^'fc') = -^E / (r-)<:(r)e^('='-'=)-[/(r) 

^ , ^ r 

~1 ^ JVi 



N 



Y.e'^'''-''>-m{r,)S^,S^6n,' J dr^-^^-^^U{r)^5,,,,{k\U{r)\k') 



where V the crystal volume, Vj the volume of the j-th unit cell, Vj the coordinate of its 
center and N the number of unit cells. We have assumed that U{r) and e*'-'^ -'=)■'' ^j-^ 
approximately constant over a unit cell and u^, (r) ~ u^if) over the range of values of 
\k' — k\ for which {k\U {r)\k') is not negligible. 

Note that usually for functions with the translation symmetry of the crystal lattice 
the scalar product is defined as 



(73) 



(*i|*2) = / dr*I(r)^'2(r) 



(with Vc the volume of the unit cell); in particular u'j^{r) and e^'^ ^'u^ are normalized 
with respect to this scalar product. 
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If the two equations (66) and (68) are identical, they have the same solutions a„ {k) and 
a{k). Thus (assuming that a,^{k) is non-zero only for the particular band n, coherently 
with our hypothesis that there is no mixing between the bands) we can write that 

(74) ^(r) = a,{k)\vk) = ^ a„(fc)^<(r) 

v.k k ^ ^ 

u^{r) ^ an(fc)^ - <(r) ^ a„{k)\k) = <(r)F(r), 



where we have assumed that u^{r) does not vary very much with k (note that the main 
fc's have to be quite close to fco = for the previous derivation to be consistent). 
We notice that if we express En{k) as 



(75) E,,{k) ^E^ + - 



m 



(with fi,u = X, z) the envelope function equation becomes 

(76) ^^(^) + + Uir)) = EF{r) 

and when the effective mass is isotropic ((1/m* = (l/m*)(5^,y) we have the weh-known 
equation 

(77) __^v2i^(r) + (£;o" + [/(r))==i?F(r). 

Luttinger and Kohn in a famous paper [40] have given an alternative derivation of the 
single-band envelope function equation, which has the advantage of being easily general- 
ized to more complicated cases. The starting equation is again the Schrodinger equation 
(i/'°-' +[/)V' = with H^^^ being the Hamiltonian of the electron in the periodic lattice 
potential and \J an additional potential which is assumed not to vary significantly over 
each unit cell. They show that the fuctions \nk) = Xnk = e^'^^'ug (r)/-\/(27r)^ (where 
Wg (r) are the Bloch lattice functions in the absence of the external potential, evaluated 
for fc = 0) are a complete orthonormal set, if considered as functions of r and k (exactly 
as the functions u^iSyV) I \/ (27r)'^, which, contrary to the Xnk-, are eigenfunctions of 
i?(o)). This means that 

(78) (Xnfe|X«'fc') = 5rm'(5(fc - fc')- 

Therefore, they can expand the wave function over the complete orthonormal set of 
functions \nk) in this way: 

(79) ^ = Yjdk'A,Ak')Xn'k', 
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and, considering this basis, they can rewrite the Schrodinger equation in the following 
form: 

(80) J2 J dk' {7jk\H^"'> + U\7j'k')An'{k') = EAn{k). 

n' 

After some calculations, they obtain that 

(81) (nfc|i/(0)|n'fc')- i?o+7^^ U„„,J(fc-fc')+ V " ° 5{k - k') = 



a—x,y,z 



where the momentum matrix elements at A; = 

(82) pr' = yj^ urHhv^W^'dr 



arc characterized by the following properties: P^™ = if the point fc = around which we 
are working is an extrcmum point of the dispersion relations, and P^*" = " — (P^" )* 
if a center of symmetry exists in the crystal. Moreover, if [/ is a "gentle" potential, with 
a very small variation over a unit cell. 



(83) {nk\U\n'k') = U{k - fc>w, 
where U (fc) is the Fourier transform of U 

(84) U{k) = j^J drer^^-Uir). 
As a consequence, eq. (80) becomes 



(85) (i5« + ^)A„(fc)+ ^ ^^^^^A„,(fc) 

a=x,y,z n' 
n 

+ j dk'U{k - fc')^n(fc') = EAr,{k). 



In order to decouple the equation corresponding to the band n from the other bands, 
the terms involving P^" , which couple the bands, have to be removed to the first order. 
Luttinger and Kohn obtain this result applying a proper canonical transformation T: 

(86) Ar,{k)=Y^ j dk'{nk\T\n'k')B„'{k'), 

n' 

which corresponds, more abstractly, to A = TB. Writing T = and applying this 
transformation to the equation (85), which can be rewritten as HA = EA, with H = 
Ha + Hb + U, we obtain {e~^ He^)B ~ EB. After some calculations, it can be proved 
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that, choosing S in such a way that Hi, + [Ha, S] = (the square brackets denoting the 
commutator), i.e. 



(87) 



hk- P"-" 6{k-k') , 
{nk\S\n'k') = { m.,iES-ES') ' " ^ " ' 
0, if ri = n' 



and neglecting the terms of order and higher and the terms which assume very smaU 
values for a "gentle" potential U, this equation becomes 



2 nie rrip 



pnn p7i n 

J2 ^<-^p J2 ik ^ 



a.^—x.y.z 



Tpn jpn'' 



Bn{k) 



n" ^n 



J 



J U{k - k')B„{k')dk' = EBn{k), 



which can be written more briefly in this form: 



(89) 



Enik)Bn{k) + ju{k- k')Bn{k')dk' = EB„{k), 



where En{k) is the dispersion relation in the absence of U{r) expanded to second order 
in k. 

Converting cq. (89) from the momentum space to the position space and defining the 
envelope function in this way 



(90) 



1 



7(2^ 

the single band envelope function equation is obtained 

(91) (-B„(-^ V) + U{r))F„ {r) = EF,-,{r), 

with En{—i'V) obtained expanding En{k) (the dispersion relation in the absence of 
C/(r)) to second order in k around fc = with non-degenerate perturbation theory and 
substituting each component of fc with the corresponding component of — i V. Being 
Fn{r) a smooth function, it has significant Fourier components only for small values of 
fe. Since for small values of k also S is small, for these components An{k) = e^Bn{k) ~ 
Bn{k) and thus, exploiting the eqs. (79) and (90), we have 



(92) 



tp^^ j dkBn{k)e 



and, noting that eq. (91) contains no interband coupling, 
(93) = F^{r)uZ{r) 
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(as already seen). If locally the external potential changes considerably within a cell, in 
such a region the equation we have derived is no longer valid, but it continues to be valid 
in regions of space sufficiently distant from it. 

Then Luttinger e Kohn adopt an analogous procedure starting from the Schrodinger 
equation written in the presence of an external magnetic field. In this way, they demon- 
strate that in such a case the envelope function satisfies an equation similar to the one 
in the absence of a magnetic field, the only difference being that the new Hamiltonian 
is obtained replacing, in the expansion of E„{k) to quadratic terms, each ka by the op- 
erator — iVc( + eAa/h (using the MKS system of units) with Aa the a-th component of 
the vector potential. Moreover in the expansion of En{k) to the second order any arising 
product of non-commuting factors has to be interpreted as the symmetrized product. 

In the case in which the extremum is at fc = feo 7^ 0, the demonstrations (both with 
and without an external magnetic field) can be repeated by just replacing MQ(r) (the 
eigenfunctions of the Hamiltonian for fc = in the absence of U{r) and of an external 
magnetic field) with 0„fcn = e^''''"^u'^g{r) (the eigenfunctions of the Hamiltonian for 
fc = fco in the absence of U{r) and of an external magnetic field). Indeed, it can be seen 
that the functions (pnK = e'^ '^''^ {e^ '^"''^ Uq (r) / y/ (27r)3) have properties analogous to those 
prevously seen for the Xnk — e'^^''^u^{r) / yJ{2TTY' considered as functions of r and /«, 
they are a complete orthonormal set of functions (such that {'~PnK.\'-Pn' n') — ^nn' 
and the momentum matrix elements computed in fco; defined as 

(94) Pr' = ^ / <o*i^k„^. - ^^^c.)uidr, 

have properties analogous to those seen in the case in which ko = 0. In this case the 
relation between the wave function and the envelope function is 

(95) ^ = F„(r)(e^^-X„(r-)) 
and the envelope function equation is 

(96) [E„ (fco - i V) + [/] F„ = , 
in the absence of magnetic field, and 



(97) 



F — EF 



in the presence of magnetic field. As before, in these expressions an expansion of £"„ 
around fco to second-order terms in — iV and in — iV -I- eA/fi, respectively, is meant. 

If there arc cxtrcma at several different values of fco within the band, we obtain 
an envelope function equation for each of them; if the solutions corresponding to the 
different fco values have different energies, the corresponding wave functions represent 
independent solutions of the Schrodinger equation; otherwise the correct wave function 
will be a linear combination of those from the different extrema associated with the same 
energy. 

When the band of interest is degenerate, Luttinger and Kohn, using a similar calcula- 
tion, arrive at a set of coupled second-order equations which correspond to the effective 
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mass equation found in the case of non-degenerate bands. In particular (assuming for 
simplicity that the degeneracy occurs at fc = 0) they assume to have, at A; = 0, r un- 
perturbed degenerate Bloch lattice functions corresponding to the same unperturbed 
energy Eg (where "unperturbed" means for fe = and in the absence of U (r) and of an 
external magnetic field) and they define them as (with j = 1, . . . , r, where r is the 
degeneracy), i.e. 



(98) = E=„c 



'J 



(notice that the (f>j's, i.e. the Uq's, can be seen as Bloch functions e^'^ '^u^ for fc = and 
thus they have to satisfy the Schrodingcr equation for fe = 0). They instead indicate as 
(with i ^ 1, . . . ,r) the unperturbed Bloch lattice functions at fe = corresponding 
to the other bands, that are not degenerate with the (j>j's. If the crystal has a center 
of symmetry, it can be proved that the momentum matrix elements between different 
0j's vanish, i.e. P^^ = 0. Luttinger and Kohn introduce the complete set of functions 
\nk) = (pnk ~ e^'''^4>n/\/ (Stt)'^ (where 0„ indicates both the (j)j's and the 0i's). Using 
this basis, they can expand the wave function in this way: 

(99) '^ = Y. I dkA^{k)4>nk 

and rewrite the Schrodinger equation as 



(100) 

thus obtaining: 



j dk'{nk\H'-"'^ +U\n'k')An'{k') = EAn{k), 



(101) (E^, + ^)A,ik)+ T.^Mk) 



a—x,y,z I 



j dk'U{k - k')Aj{k') = EAj{k) 



(writing only the equations corresponding to the degenerate states j). In order to de- 
couple the equations corresponding to the states j from those of the states a proper 
canonical transformation A = TB = e^ B is again applied, with 

r hk-p^'^'sik-k') , . 

(102) (nfe|5Kfe')= - nMES~ES') ' " ^ " ^ ^ ''l' 

[ 0, if n and n' G [1, r]. 



In this way Luttinger and Kohn obtain, to second-order terms in fc, the following set of 
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equations for the r degenerate states: 

(103) j2\^oSn'+ E {Djj^'k^h)] Brik) 

j' = l y a,0=x,y,z j 

+ y" U{k - k')Bj{k')dk' = EBj{k), 



with 



Therefore, introducing again the envelope functions 

(105) F,{r) = ^=L= J e^''-B,{k)dk, 

Luttinger and Kohn arrive at the conclusion that the r envelope functions Fj{r) cor- 
responding to the originally degenerate energy bands satisfy the r coupled differential 
equations 

(106) I ^0%' + E + Uir)S,,, J F,,{r) = EF,{r) 

j' = l y a,p=x,y,z j 

(if the energy zero is set at the term E^^jy disappears). 

Analogously to what happens in the non-degenerate case, for small values of fc, 
A„(fc) ~ Bn{k) and thus 

(107) dkB^{k)e^^-^^ = Y,FrXr)K{r) c^j2F,{r)<j,,{r), 



since in eq. (106) no coupling remains between the states j and the states i. The num- 
bers Dj^, play the same role in the case of degenerate bands as ?i^/(2m*^) for a non- 
degenerate band. 

As before, in the presence of a magnetic field the components of — iV appearing in 
the envelope function equations will be replaced with the corresponding components of 
-iV + eA/h. 

In the presence of spin-orbit coupling, Luttinger and Kohn adopt the same treatment, 
considering the spin-orbit contribution as part of the unperturbed Hamiltonian (therefore 
the total unperturbed Hamiltonian will be 77^°^ -I- Hso) and assuming the Bloch lattice 
functions and the corresponding energies for fc = of H^^^ + Hso as known quantities. 
Thus the Ug are replaced with the Mq (the spinorial Bloch lattice functions for fc = in 
the presence of spin-orbit interaction), i?„(fc) by En{k) (the dispersion relations in the 
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presence of spin-orbit interaction) and the " by 

(108) (7r„)„„, = (u^l {-ihVa, + -r^{(T x VV)^ K'), 

where the extra term arises from the fact that the spin-orbit couphng contains the dif- 
ferential operator p. When we treat energy bands which are degenerate in the absence 
of spin-orbit interaction, we have to remember that (as seen previously) the spin-orbit 
coupling can hft, at least partially, the degeneracy. In such a case, we have to consider 
that the validity of the adopted theory rests on the assumption that the interband sep- 
arations are large compared with the energies involved in the solution of the envelope 
function equation. Thus we have to evaluate if the external potential \J or the magnetic 
field are sufficiently small to produce no appreciable mixing of the bands, the degeneracy 
of which has been lifted by the spin-orbit coupling. If they are sufficiently small, we can 
obtain a different set of coupled envelope function equations for each set of bands that 
have remained degenerate; otherwise we will have to deal with the full set of coupled 
equations for all the bands that are degenerate in the absence of spin-orbit. 
We can introduce a matrix Z?, the elements of which are 

(109) D,,. ^Y.^t^-^?- 

If in these matrix elements we replace each component of the vector k with the corre- 
sponding component of the operator — i V -I- eA/h, we obtain the terms which appear 
in the envelope function coupled equations. In particular, the envelope function coupled 
equations written in the absence of an external perturbation read (if we set the energy 
zero at E^) 

r 

(110) J2 E(^«^(-*^")(-*^^))^:'' W = ^^^' W- 
j'=l a,P 

If we convert them from the position representation to the momentum representation, 
we obtain 

r 

(111) Y.^D';PKk,)B,,{k) = EB,{k) ^ 

j' = l a,l3 
r 

Djj^Bj,{k) = EBj{k) ^ DB = EB, 

i'=i 

from which it is evident that the dispersion relations E(k) near the extremum can be 
obtained by finding the eigenvalues of the matrix D. We notice that this clearly corre- 
sponds to what happens in the case of non-degeneracy, in which (as we have seen) the 
envelope function equation contains £'„(— i V) (the dispersion relation in the absence of 
external potential energy or magnetic field, where each component of k is replaced with 
the corresponding component of — i V). 

In order to determine the number of independent parameters which appear in the 
matrix D, the symmetry properties of the considered lattices are exploited. 
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Fig. 3. - Heterojunction between two semiconductors A and B. 



In [41] Luttinger proposes a different way to obtain an explicit expression for D, based 
only on symmetry arguments. He writes this matrix for diamond-type semiconductors 
using group theory, in particular considering that the Hamihonian D should be invariant 
under the operations of the cubic group (so that the Hamiltonian will give us results 
which transform correctly with respect to the transformations of the cubic group, which 
is the symmetry group of k) and thus writing _D as a linear combination of the invariants 
obtained combining angular momentum matrices and components of k. The elements 
of such a matrix are polynomials in the components of fc, at most of the second order, 
and involve parameters characteristic of the materials, which have been experimentally 
found and are available for most common semiconductors [79]. For example, in the case 
of the 4x4 matrix D corresponding to the light-hole and heavy-hole bands (the extra 
factor of 2 coming from spin) they are 71, 72, 73, k (which is useful in the presence of an 
external magnetic field) and q (which approaches zero as the spin-orbit coupling does). 

Bir and Pikus [14] have shown that in uniformly strained semiconductors, such that 
the periodicity of the structure is preserved, the strain introduces in the dispersion rela- 
tion of non-degenerate bands an extra term of the kind 

(112) adexx + ^yy + f^zz) 

and in the Hamiltonian of degenerate bands additional terms of the form 



where a, 13 = x,y, z and Sap is the generic component of the strain matrix. 

Bastard [42-44] uses the envelope function method to study heterostructures, for 
example made up of two materials A and B (fig. 3). In particular, he assumes that the 
two materials are perfectly lattice-matched and crystallize with the same crystallographic 
structure, so that the functions Ugir) in the two materials can be considered identical. 
With this hypothesis, if in each material the wave functions are written as 



(113) 




(114) 
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it is evident that, since the Uq are Hnearly independent and the wave function has to 
be continuous at the interface, also the envelope functions have to be continuous at the 
interface. For the derivative of the envelope functions, Bastard finds, enforcing the con- 
tinuity of the probability current at the interface, a general condition [80] , which, in the 
simple case of two materials that are both characterized by non-degenerate parabolic and 
isotropic bands but with different effective masses "i*^^ and 'm*^Qy reduces to enforcing 
the continuity of 



(where we have assumed the z axis orthogonal to the interface). This can be easily 
obtained enforcing in this case the continuity of the z component of the probability 
current density, which is equal to 



and noting that the continuity of the envelope function has already been enforced. As to 
the asymptotic behavior of the envelope functions far from the interface, it depends on 
the heterostructure under consideration. For example, for superlattices the z-dependent 
part of the envelope function will be a Bloch wave, due to the periodicity of the structure 
in that direction, while for the bound states of a quantum well it should tend to zero for 
large z. Thus the envelope functions in the overall structure can be found solving the 
envelope function equations in the different materials, knowing the asymptotic behavior 
far from the interface and enforcing the correct boundary conditions at the interface. 
Bastard has also made an extensive analysis of the applications of this method [81]. 

Also M. Altarelli has given important contributions to the development of the envelope 
function method [82] and to its applications to the study of hetcrostructures [45-47]. 

M. G. Burt [48-52] has pointed out the errors deriving from the assumption, normally 
made in the application of the envelope function method to hetcrostructures, that the 
Mg(r) in the two materials are the same and from the boundary condition enforced on 
the derivative of the envelope function at the interface. In a series of interesting and 
detailed articles he has developed an alternative envelope function theory expanding 
the wave function in the overall structure on the same periodic basis functions Un{r) 
throughout, even though they are not necessarily eigenstates of the constituent crystals, 
without making any hypothesis about the real eigenstates Uq (r) 



The envelope functions F„{r) univocally defined in this way and all their derivatives 
are certainly continuous everywhere, including at the interface. Using this approach, 
he has first derived exact envelope function equations, then, for local potentials and 
slowly varying envelope functions (but without any assumption on the rate of variation 
of the composition), he has formulated approximate envelope function equations, and 
finally, with the assumption of the dominance of one envelope function, he has arrived 
at an effective- mass equation that includes also the effect of the differences in the (r) 



(116) 




(117) 




n 




between the two materials. At each step the associated approximations are accurately 
described, so that it is possible to estimate the error. 

A more detailed description of the applications of the k ■ p method to materials with 
a diamond, zincblende and wurtzite lattice, both in the periodic and in the non-periodic 
case, can be found (besides in the other books and in the original publications reported 
in the list of references of this review) in the recent book by L. C. Lew Yan Voon and 
M. Willatzen [4]. 



4. Application of the k ■ p method to graphene 

In the last years the k ■ p method, and in particular the formulation (described in 
the last section) based on the envelope functions, has been successfully applied to the 
analysis of the electronic properties of graphene and graphene-related stuctures, such as 
carbon nanotubes and graphene nanoribbons. 

In this section we will begin the description of this particular application deriving the 
k ■ p relations for a simple sheet of graphene. 

A graphene sheet is a hexagonal lattice of carbon atoms. In fig. 4(a) we show its 
structure in the real space and, in particular, its unit cell as a dashed rhombus, containing 
two inequivalent carbon atoms A and B, while in fig. 4(b) we show the lattice in the 
reciprocal space with the Brillouin zone as a shaded hexagon. The lattice unit vectors 
are ai and a2 in the real space, and bi and 62 in the reciprocal space. If we define 
a = \ai\ = \a2\ = ac-cV^ (with ac-c the distance between nearest-neighbor carbon 
atoms), the coordinates of these vectors in the right-hand reference frame S' = [x' , y' , z') 
are (observe that we have taken x along the vector ai -1-02) 



(118) 



ai 



V3 



0.2 



V3 



bi 



2n 

\/3a 
27r 

a 




2tt 

\/3a 
27r 







(following the conventions used by R. Saito, G. Dresselhaus and M. S. Dresselhaus [83]), 
which (being 61 — 27r(a2 x z')/(ai • (02 x z')) and 62 = 2n{z' x ai)/{ai ■ (02 x z'))) fulfill 
the well-know relation • bj = 2'KSij between lattice unit vectors in the real space and 
in the reciprocal space. Note that the letter written under the symbol "=" indicates the 
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Fig. 5. - The energy dispersion relations of graphene inside its hexagonal Brillouin zone. 



adopted reference frame. The most relevant graphene dispersion relations for transport 
and other solid-state properties are the two 7r-bands (an upper anti-bonding band and a 
lower bonding band), which are degenerate at the points (considering the point T at the 
center of the hexagonal Brillouin zone of graphene as the origin of the reciprocal space) 



and obviously at their equivalents in the reciprocal space (as we can see from fig. 5, 
which has been obtained with a nearest-neighbor tight-binding approach limited to the 
2p^ atomic orbitals, with nonzero nearest- neighbor overlap integral). 

Thus we can use the k-p method to find the dispersion relations of graphene near these 
extrema points (called Dirac points), following T. Ando's approach [54-56]. However, in 
our description we will continue to use the conventions of ref. [83] and we will consider 
the pair (119) of Dirac points (which will simplify the treatment of zigzag and armchair 
graphene nanoribbons in the last section of this review). Other articles where a, k ■ p 
treatment of graphene is introduced are refs. [84-89]. 

We start using a simple tight-binding model, in which we use as basis functions the 2p^ 
orbitals of all the carbon atoms of the graphene sheet, which are the orbitals leading to 
the TT-bonds and thus to the above-mentioned two 7r-bands. The generic eigenfunction in 
the material can be expressed [55,56] as a linear combination (with coefficients iI}a{Ra) 
and 'iPb{Rb)) of these atomic orbitals </3(r — Ra) and <^{r — Rb) (centered on atoms of 
type A and B, respectively) 



(119) 



K = 




(120) 



V'(r) = ^ i,A{RA)^p{r ~Ra)^Y^ ^JBiRsMr ~ Rb) 




where the first (second) sum spans over all the positions of the atoms of type A (B) in 
the lattice. 
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(121) 

we have that 



HM=E[ 



(122) 



and thus (using j and j' to indicate the type of the atoms and n and m to specify the 
particular atoms) 



(123) E 



\j=A,B Rj„ 



H 



( H 51 Tpj{Rjn)v{r ~ Rjn) 51 51 "i^y {Rj',„)ipir - Rj.rn) ) 
\j=A,B Rj„ j'=A,BR^,^ I 

J2 r,{RM'{Rfra){v{r - R,n)\H\^{r - R,,ra)) 

jJ'=A,B Rj„ R-,^ ^ 

X! X! i^jiRjn)^j'iRj'm){v{r - Rjn)\'pir - Rfm)) 
j,j'=A,BRj„ R.,^ 

E E E V'*(i2,n)V','(i2,'™)/iR,„,R,,„ 
j,j'=A,BRj„ 

E E E ^i(-Rj»)V'j'(-Rj'm)sH,„,R^.,,„ ' 
jJ'=A,BR,„ R.,^ 

where we have introduced the transfer integrals hnj^^R.,^ and the overlap integrals 
SRj„,Rj'm between atomic orbitals. Now we can minimize E (to obtain the actual physical 
state) enforcing (for each coefficient, and thus for each atom) 



(124) 



E E V'j'(iij'm)/lR,„,R^,„ 
dE _ j'=A,BRj,^ 

j,j'=A.BRj„ R-,^ 

E E E rj{R,n)^r{Rrm)hR,,..R, 

j.j'=A.BRj„ R-,^ 



E E E V'j(-Rj«)^/'j'(i2j'm)sR,„,R^., 
j.j'=A,BRj„ 
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Multiplying both members by the denominator of eq. (123) and rearranging, we find: 

(125) T.'f'riRj'r.^)hR,„,R,,^ = 

3'=A,BR.,^ 

E E E V'l(-Rjn)V'i'(^j'™)'lH,„,H,,,„ 
j,j'=A.B Rj„ Rji,„ \ ^ . ^ 

.r-^ 2^ 2^ i^]'{Rj'ni)sR^„,R ,^ 

'm)SR-^.R.,^ ji =A,B R-, ^ 

3,3'=A,BRj„ R.,^ 

and recognizing that the fraction in the right-hand side is the expression of i?, we have 

(126) Y.'f'riRrrn)hR^„M,,^=E J2 J2^r(^rm)sR,„,R^,„^. 

]'=A,BR^,^ ]'=A,BR.^ 



Let us expand this result for the coefficients (and thus for the atoms) with j — A and 
for those with j = B 



(127) 



E '^A{RA7n)hR^„,R^^ + Y '>pB{RBm)hR^^.Rs^ = 

RAth Rsm 

Am)sRj^^.RAm+'Yj ''l^B{RBm)sRA^.RB 

\Ra^ Rsm 

Y '>pA{RA7n)hRs^,R^^ + Y ^B{RBm)hRs^.Rs„^ 

RAtu Rsm 

E Y^ i^A{RA7n)sRs„,RAm + Y. '^B{RB7n)sRs„,RE 
Vka™ Rb^ 



We consider non-negligible only the integrals between each atom and itself and between 
each atom and its nearest neighbors (which are the nearest three B atoms for an A 
atom, while they are the nearest three A atoms for a B atom). Therefore, if (in order to 
simplify the notation) we rename Rau as Ra and Rbu as Rb and wc use the index I to 
indicate the nearest three atoms, we can rewrite these equations in the following way: 



(128) 



ijA{RA)hR^.R^ +Y'^B{RB,)hR^,Rs^ = 

1=1 

E (^AiRA)sR^,R^ +^V'B(-RB,)sfl^,flB, j 

3 

Y'^A{RAi)hRs.R^^ +i}B{RB)hRg,Rs = 

1=1 

^ {^Z^a{RAi)sRs,Ra, +i^B{RB)sRs.R^ 
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In particular, we consider 



(129) 



-70 




u{Rjn) if Rj 



1 if Rjn — Rj' 

if Rjn Rj' 



Rj'jni 



if Rjn 7^ Rj' 

otherwise, 



and Rjn and Rj'm are 
nearest neighbors. 



Here 70 is the modulus of the nearest-neighbor transfer integral. Instead tR^^ is the 
onsite energy, that wc take as zero of the energy in the absence of an external (i.e. not 
due to the periodic structure of the lattice) potential energy; if the external potential 
energy is not zero, we have to consider the term u{Rjn), which represents the value of 
this external potential energy in the position Rjn- 

Note that the reason for the values of the overlap integrals reported in eq. (129) is 
that we consider atomic orbitals orthonormalizcd using the Lowdin procedure [90-92]. 

Thus the tight-binding relations become 



(130) 



-joY^i^siRB,) = {E-u{Ra))^a{Ra), 

1=1 
3 

'10^^a{RA:) = {E-u{RB))i^B{RB). 
1=1 



If we introduce the vectors (fig. 4(a)) 



(131) 



T2 



1 
2 

2 




T3 



1 
2 

V3 
2 




(with respect to the frame S' = {x' ,y' , z')), we can write the positions of the nearest- 
neighbor atoms in this way: 



(132) 



Rbi = Ra — Ti , 

Rb2 = Ra — T2 , 

Rbs = Ra — T37 

Rai ~ Rb + Tl, 

Ra2 ~ Rb + '''2, 

Ras = Rb + '''3, 
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and thus we can rewrite the tight-binding relations in the following form: 



(133) 



-70 J2 ^b{Ra - Ti) = {E- u{Ra)) ^a{Ra), 
1=1 

3 

-loY^i'AiRB + Ti) = {E-u{Rb))4'b{Rb). 

1=1 



Now let us consider what happens near the points K and K' . 
Let us assume that we can write 



(134) 



^B(flB) = »e'^'e'^-^-i^f (Rb) + e'^'-«-F|^'(fiB) 



(the angle 6' will be properly chosen later). If fc is the wave vector of ipA and tpB, 
the functions and have a wave vector k = k — K and thus are slowly- varying 
functions (with small k) near the point K; analogously the functions F^ and Fg have 
a wave vector k = k K' and thus are slowly-varying functions (with small k) near the 
point K' (note that in the overall review we use fc for the total wave vector and n for 
the wave vector measured from the reference extremum point). 

Incidentally, with these assumptions, if we define = 1, 
and = 1, we have that 



I e 



10' 



(135) ij{r)= ^ ^V^(i20'/'(r--«0 = 

i=A,B Ri 

EE E «^e'^-^-^;^^(ii.)^(r--i2.) = 

i=A,B Ri Kj=K,K' 

EE E «^e'^-«-F/^^(r)^(r.-l?,) = 



i=A,B Ri Kj=K,K' 



i=A,B Kj=K,K' 



Ri 



i=A,B Ki=K,K' 



K K ■ K ■ 

where we have substituted F^ ' (r) to F^' (Ri) using the fact that F^ ^ is a slowly varying 
function of r near Kj, while the atomic orbital (p has significant values only near the 
corresponding atom. The quantity between square brackets (that we have called here 
, ) is periodic with the periodicity of the lattice, since, if a£ is a lattice unit vector 
(and thus, if Ri is the position of a lattice point, also i?^ = Ri — ai is the position of 
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a lattice point), then 

(136) iii^^(r + a£) =«f^-^^((r + a,)-il,)e-^^^-«'^+"^)-^-) = 

Ri 

af^- ^^(r - [R, - a,))e-^^^-('^-(«'-"^» = 
af^$:^(r-i2°)e-^-(-«")=.iW^(r). 



Therefore, since . has the lattice periodicity and Kj is an extremum point (different 

from 0) of the dispersion relations, from the relation between tp{r), '!i^ .(r) and ^(r) 

we conclude that the 4 functions F^ ' can be seen as the electron envelope functions 
corresponding to the 2 extremum points Kj where the 2 considered bands of graphene 
are degenerate (see eq. (95), the related discussion, and eq. (107)). 

Let us point out that this whole procedure does not need a particular choice of scalar 
product and of normalization: these have just to be chosen coherently with each other. 
However, one could find desirable to normalize the periodic function . according to 
the scalar product defined in (73), as is generally done in the envelope function theory. 
Following this particular criterion, one should have (if i7o is the area of a graphene unit 
cell, while Q, is the area of the overall graphene sheet) 

(137) l = (S^.(r)|^^^^.(r)) = -l / |^k,Wpdr=l/ MP^r = 

1 / ( ^ ^(r - R^) e-^H-R.)) * ( ^(r - H,) e-^H-^D) rf^ ^ 



Ri 



R' 



iKy(R[-R^) 



^l^\,ir-R.)\^dr 

+ A E V*ir-R,Mr~R!;)dr 

" R.M 
Rz^R'i 

1 / J2\^ir-R,)\'dr^^J2 I W{r^R,)\^dr 



Here we have exploited the following properties of the involved functions. First of all, 
integrating a function with the lattice periodicity over the whole graphene sheet and 
dividing the result by its area is equivalent to integrating it over the lattice unit cell and 
dividing by the corresponding area. Moreover, each atomic orbital ^p (orthonormalized 
using the Lowdin procedure) has a non-zero overlap only with itself. Finally, since each 
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atomic orbital has significative values only near the corresponding atom, the integral of 
the square modulus over the whole graphene sheet is nearly the same for all the considered 
atomic orbitals, and thus the sum of all the integrals is approximately equal to a single 
integral multiplied by the number ri/Oo of orbitals. 

Therefore, adopting this particular normalization for . , the atomic orbital ip should 
be normalized in such a way that 

(138) / \^{r ~ R,)\^dr = 1 ^ / |^(r - R^)\^dr = f^o, 
"0 jvi Jn 

and thus we should consider atomic orbitals v^Hq times greater than those deriving from 
the usual normalization over the whole graphene sheet. 
The corresponding scalar product 

(139) (^i|^2) = ^ / ^t{r)Mr)dr 

"0 Jn 

should be used in all the calculations involving atomic orbitals. 

If we introduce the assumptions (134) into the tight-binding equations (133), we 
obtain 



(140) 



1=1 



{E-u[Rb)) 



le e 



'JK-Rb TpK 



B) 



JK' Rb TpK' 



F§ {Re) 



[e*^-(^"+"')<(i2i3 +Ti)-i e'''e''^'<''^+^''>Ff{RB + n) 



'70 > le 

1=1 



It is useful to introduce [55,56] a smoothing function g{r), i.e. a real function which 
varies smoothly around the point around which it is centered, has non-negligible values 
only in a range of a few lattice constants around the center, and then decays rapidly for 
larger distances. This function (point-symmetric around its center) is chosen in such a 
way as to satisfy the conditions 

(141) ^5(r-fi^)==^g(r-JlB) = l 

Ra Rb 



and 
(142) 



/ drg{r-RA)= / dr gir - Rb) = 
Jn Jn 



(where Oq = ■\/3a^/2 is the area of a graphene unit cell, while Vt is the area of the overall 
graphene sheet); moreover it has to satisfy the relations 

(143) - RA)e'^'^'~'^^-''^ ^J29ir~ fiB)e'(^'-^)-«« ~ 0. 

Ra Rb 
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1 — , 



g(x.y) 0-5- 




Fig. 6. - A candidate smoothing function g{r). 



Due to its locality, when this function is multiplied by a generic smooth function f{r) 
(such as the envelope functions F we have defined), we clearly have that 

(144) /(r).g(r - R) f{R)g{r - R) 

(for positions r for which g{r — R) is not negligible, the smooth function /(r) is ap- 
proximately equal to f{R), while for positions r, further away from R, for which f{r) 
significantly differs from f{R), the function glr^-R) is null). In fig. 6 we show a possible 
smoothing function g{r), which approximately satisfies all the previous relations (^). 

If we multiply the first of the tight-binding equations (140) by g{r — RA)e~^^'^-'^ and 
we sum it over Ra, we find 

(145) EY,9{r~RA)F^{RA) 

Ra 

-E^^<^'Y.g{r-RA)e^^^'-^y^-Ff{RA) 
Ra 

- 9{r - Ra)u{Ra)F^ {Ra) 
Ra 

+ie''' Y^g{r ~ RA)e''^'''-''y''-u{RA)Ff [Ra) = 
Ra 

-70 I e^'' E E 9{r ~ Ra)F^{Ra - rz) 

1=1 Ra 
(=1 Ra 



(^) In detail, we have represented the function defined as 106.5307 exp ( — ^ _ ^ | ^ | /^('J^ 355 nm) ) ) 
|r| < 0.355 nm, and for \r\ > 0.355 nm, but better approximations for the smoothing function 
g{r) can be found. 
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exploiting the property (144) it becomes 



(146) E 



F^{r)-Eie'' 



LRa 

Y,9ir~ RaHRa) 



LRa 



Ffir) 



Ra 

5^g(r-H^)e*(^'-^)-«-u(i2. 
Ra 

^e^^Xe-^- \j2g{r-RA) 



Ff{r) 



-7o I e > e 

1=1 

3 

-iK'-i-i 



LRa 



F^{r-Ti) 



-10 J2' 

1=1 



Y,g{r~RA)e^'^^'-^y^- 
Ra 



F^{r-Ti). 



For the quantities in the square brackets, we can use the properties (141) and (143), 
together with the definitions 

(147) UA(r) = ^5(r-i?A)^(fiA), u'j,{r) ^Y^di^ ^ ^^y'^'^' 
Ra Ra 



obtaining 
(148) 



E F^ (r) - UA{r) (r) + » e"^' u'^r) F^ (r) 

3 

-Jo^e'''y^e-^^■-'F^ir-Tl). 



?K', 



1 = 1 



Expanding the smooth quantity Fg{r ~ ti) to the first order in r;, we have that 

3 3 



(149) ^e-'^^-'i^f(r-r0^5], 



-iK-Ti 



1=1 



1=1 



F^{r)-[n.^]F^ir) 



\i=i 



-iK-T, \ ipK 



F^ir)- 



Let us now calculate the value of the sums which appear in the previous expression 

3 

(150) ^e^'-f^-^' = l + g-^T +e*T ^0; 



1=1 

3 



1=1 



Tl 



d 



= 1- 



a / d 



19 V3 a ' 



_j2x a 

' V§ I 2 9a;' 2 dy' 



Id d 



a 

7! 



1 _-2^ 1 2^ 

-1 + -e ' 3 + -e* 3 

2 2 



'^^ ^ ^\2dx' ' 2 9y' 



dy' 
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1 _,2lL 1 3 , \/3 _,,27r \/3 j2n_ ,3 , , 

Since — l + -e 3 + -e 3 = and e 3 ^ g 3 = i-, we have that 

2 2 2 2 2 2' 

(151) te-'^-'(..l]^~A.U^-^' 



1=1 



dr J V32 \dx' dy' 



where we have defined = — iV and thus 



(152) kx' = and 



d , . . d 



" dx' ' ^ dy'' 

Substituting these results, eq. (148) becomes 

(153) EF^ir) - UA{r) F^{r) + ^ e'"' u'j,{r) Ff{r) 

-loie'^ i i^a{kx' ~ iKy')Fg {r)\ ^ 



^Joae^'^' {kx' - iky,)Fg{r) = -/{k^ - iky)Fg{r), 

where we have passed from the original reference frame S' — {x , y , z') to a new frame 
S = {x, y, z), rotated, in the plane {x , y), around the origin by an angle 0' (positive in 
the counterclockwise direction) with respect to the original one (fig. 7) and we have used 
the fact that 

(154) e*^ {kx' — iky') = {cos0' + isin0'){kx' — iky') = 
{cos0'kx' + sin0' ky') — i{cos0' kyi — sin0' kx') — kx — iky 

(due to the relations between old and new coordinates), with 

d d 

(155) kx = —i-^ and k„ = —i-z—. 

^ ' dx ^ dy 

Indeed, it is a well-known result that, for a rotation by 0' of the reference frame, the 
relations between the new and the old coordinates are x ~ x' cos 0' + y' sin 0' and y = 
y' cos 0' — x' sin 0' . Therefore we have that 

(156) ^■F(a;,y) ^ dF{x,y) dx ^ dF{x,y) dy ^ dF{x,y) dF{x,y) 

dx' dx dx' dy dx' dx dy 

and that 

(157) g-F(a:,y) ^ dF{x,y) dx ^ dF{x,y) dy ^ dF{x,y) ^ dF{x,y) 

dy' dx dy' dy dy' dx dy 
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Fig. 7. - The reference frames used in the calculations (Ch and 9 will be used for carbon 
nanotubes in the next section: this figure corresponds to a (4, 2) nanotube). 



As a consequence, we have that 



(158) {cosO'k^' + sin6'ky,)F{x,y) = cos6'' ^- 



. dF{x,y) 
dx' 



: dFix,yy 



dx 



dy 



dF(x,y) . 9 dF{x,y) . 

^ ^^sin^e^'H sin 6*' cos 



dx 
. dF{x,y) 
dx 



(cos^ e' + sin^ e') 



dy 

2n'\^ ,dF{x,y) 



dx 



KxF{x,y) 



and that 



(159) (cos6I'kj,. - sine' k^,)F{x,y) = cosO' {^i ^^'^f^ ^ - sing' ^- 



, dF{x,yy 
dx' 



9Fix,y) , dF{x,y) 



dx 
dF{x,y) 



sma COS! 



„2 01 



dy 



dx 

. dF{x,y) 
dy 



COS t/ sm 



r + ^^:p!l sin' 0' 
dy 



(cos^ 9' + sin" e') = = k.,F{x, y), 

dy 



from which we obtain eq. (154). 

9' is the angle, taken counterclockwise, from the vector ai + 02 to the axis x of the 
new frame. We have also defined the quantity 7 = ('\/3/2)7oa. 
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Note that in the new reference frame S = {x, y, z) 



(160) ai =E 



2-K 



V3cos6'' + sm6'' 
008 6*' - VSsine*' 



cos 6*' + V3sin6'' 



0.2 



s 2 



3a 



cos a — sni ( 


-sine' ' 
-cose' 




2ti 



, 47r 
K' = — 

s 3a 



%/3cos6l' - sin6l' 
-cos 61' - V3sin6'' 


cose' sine' 
-\/3 cos 6*' - sin 61' 


sin 61' 
cos e' 




Analogously, if wc multiply the second of the tight-binding equations (140) by g{r — 
RB){—ie~^^ e~'^^'^B^ ^^t^ ^y^jn it over Rb, using again the properties (144), (141) 
and (143), together with the definitions 

(161) UB{r)^Y.9i^'^B)u{RB), u'B{r)=Y,g{r^RB)e'^'^'-'^^-^-u{RB), 

Rb Rb 

we obtain [93] 

(162) EF^{r)-UB{r)F^{r)+ze-^o' u'Bir) Pf (r ) = 

3 

1=1 

Expanding the smooth quantity F^{r + t;) to the first order in t;, we have that 

3 3 



(163) ^e'^^-'Fi^(r + r,)^5]e*^--' 



(=1 



1=1 

3 



dr 



.1=1 



Since 



(164) J2e'^-^^=0 and ^ e'^-' • ^ = -*^a(/i,, + 



d 



.^/3 



dr 

1=1 1=1 
eq. (162) becomes 

(165) EF^{r)-UB{r)F^{r)+ie-^''u'B{r)Ff{r) 
-ie' 



70 I e 
a/3 



-i^a{k^, + iky,) \ F^{r) 



7oae (k^- + iky,)F2 (r) = 7(K:r + iKj,)Ff (r) 
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where we have made use of the relation 



(166) 



e {kx' + iiiyi) = {cos9' — i sm6'){fix' + iHyi) = 
(cos^'k-e' + s'mO'fiyi) + i{cos9'ky/ — siii9'k.j;>) ^ Kx + 



Instead, if we multiply the first of the tight-binding equations (140) by g{r — Ra) x 
{i e~*^ e~^^ ) and wc sum it over Ra, wc obtain (exploiting the properties (144), (141) 
and (143)) [93] 



(167) 



EFf{r)^^e-^'>'u'A*ir)F^ir)-UA{r)Ff{r) = 



-id' J 



?K' , 



1=1 

Expanding the smooth quantity Fg {r — Ti) to the first order in T;, we have that 

3 3 

(168) Yl e-'^''"'FB' [r - r,) ^ ^ e''^'-^' 



1=1 



\i=i 



-iK'-Ti pK' 



(r) 



1=1 

3 



.1=1 



d_ 
dr 



with 



(169) =0 E V'-Q^)^ -i^a{kx' + iky,) 



1=1 



Therefore eq. (167) becomes 



(170) 



EFf{r) ~ie-^''-u'/{r)F^{r) ~ UA{r)F^' (r) 
-70^6-*'^' (i^a{kx' + iky,)Ff {rU = 



1=1 



-iS',' *f„\ ipK 



F^'ir), 



.V3 



K'l 



—-foae~'^\kx' +iky,)F^'{r) = j{kx + iky)F§-' (r), 

where we have exploited the relation (166). 

Finally, if we multiply the second of the tight-binding equations (140) by g{r — 
Rb) X e^'-*^ and we sum it over Rb, we obtain (using the properties (144), (141) 
and (143)) [93] 



(171) 



EF^{r)-z e'O u'b* [r) F§ [r) ^ ub (r) F^ (r ) 



70 « 



1=1 
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Expanding the smooth quantity {r + ti) to the first order in t; , we have that 



(172) e''^''"'Ff (r + r,) ^ ^ e'"^''^' 



1=1 



\i=i 



1=1 

■ 3 



' dr 



Ff{r). 



Since 
(173) 



^e*^'--'=0 and ^ e'^''^' ( • ^ ) = -z^a(K,- - m^O 



1=1 



1=1 



eq. (171) becomes 
(174) 



E (r) - ^ e^' i4'(r) F^ (r) usir) F§ (r) = 



7K' I 



70 I e 
V3 



-i—a{k^, - iky,) ) Ff'(r) 



700 6'" (kj;' - iky')F^ (r) = 7(k^ - iky)F^ (r), 

where the relation (154) has been used. 

In this way, we have obtained the four equations (153), (165), (170) and (174), that 
we can summarize 



?K' 



K' I 



{ UA{r) F^{r) + 7(^. - iky)F^ {r) -i^' u'^r) Ff{r) =EF^{r), 



?K', 



7(ac, + tky)F^ (r) + UB{r) F^ (r) - z e'^'' u'^ir) F^' (r) = E F^ (r), 
te-^''u'/ir)F^{r)+UA{r)Ff{r)+^{k,+tky)Ff{r) = EFf{r), 



?K' , 



(175) <^ 



L ^e^'' u's*{r) F^ (r) + j{k, ^ky)Ff{r) + usir) F^' [r] ^ E F§' (r) 
and write in matrix form 



7K' , 



(176) 



UA{r) 
j{kx + iky) 
^e-*«' u'/{r) 




E 



Ff{r) 
F^{r) 
Ff{r) 
Ff{r) 



j{kx-iky) -ie''^ u'^{r) 







usir) -ie-'" u'^{r) 

UA{r) l{kx + iky) 

ie"^' u'g*{r) -yikx-iky) usir) 



(r) 
Ff (r) 
Ff{r) 
Ff{r) 



which is the fc • p equation of graphene. 
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Incidentally, if we repeat all the previous calculations considering the following differ- 
ent pair of reference Dirac points: 



(177) 



K 



- 27r " 




- 27r ■ 


\/3a 






27r 


, K' = 


27r 


3a 




3a 











(equivalent, in the reciprocal space, to the pair (119) of Dirac points), wc have to replace 
(134) with 



(178) 



where rj = 7r/6 + 0' , and we obtain (instead of cq. (176)) 



(179) 



E 



UA{r) 
^{k.^ + iky) 



(r) 
Ff (r) 
Ff (r) 
■(r) 



-f{kx~iky) e*''-u^(r) 

UB{r) 

UA{r) 

^e*^ e"' u'g*{r) jik^ ~ iky) 





j{kx + iky) 
UB{r) 



F^{r) 
F^{r) 
Ff{r) 
Ff (r) 



as found by Ando [55,56]. 

Summarizing, wc have that the overall wave function is given by (sec (120)) 



(180) 
with (see (134)) 
(181) 



^(r) =J2^A{RA)vir - Ra) 
Ra 



'ijjBiRB)^p{r " Rb), 



i^Air) 
il>B{r) 



^F^{r)~ie^''e^'''-^F^\r), 
'e^^-'-Fj^(r) + e'^'-Ff' (r). 



where the envelope functions F satisfy cq. (176). 

We can treat two limiting cases for the external potential, depending on its range [94, 95] . 

If the potential range is much smaller than the lattice constant (short-range case), we 
can consider the external potential as different from zero only on one carbon atom. 

If it is non-zero only on an atom of type A (in position Raq), i-e. u{Rao) ^ 0, 
u{Ra) — for Ra 7^ Rao a-nd u{Rb) = for every Rb, recalling cq. (147) and (161), 
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we have that 

(182) UA{r) = J^air - RaHRa) = g{r - RaoHRao), 

Ra 

u'Air) = E - RA)e^'^'^'-'^^-''-u{RA) = 

Ra 

UB{r) = E 5(r - Rb)u{Rb) = 0, 

Rb 

Mr) - E - RB)e^^''' u{Rb) - 0. 

Rb 

Instead; if it is nonzero only on an atom of type B (in position Rbo), i.e. u{Rbo) 7^ 0, 
u{Rb) ~ for Rb 7^ Rbq a-nd u{Ra) = for every Ra, we have that 

(183) UA{r) = J29ir-RA)uiRA)=0, 

Ra 

u'Air) = Y.9ir- RA)e'^'''-''^-''MRA) = 0, 
Ra 

'^B{r) ^^g{r- Rb)u{Rb) g{r - Rbo)u{Rbo), 

Rb 

u'Bir) ^Y.9ir~ RB)e'^''' u{Rb) = 

Rb 

g{r - i2so)e'(^'-^>^«oy(HBj = tiB(r)e'(^'-^>««o . 

If instead the potential range is much larger than the lattice constant (long-range case), 
using cq. (141), (143) and (144), we have that 



(184) UA{r) = Y,9{r - RA)u{RA)^Y.3i^ " RA)u{r) 



Ra 



Ra 



Y.9ir~R^ 



Ra 



u{r) — u{r), 



<{r) = Y,9{r- RA)e'^'''-''>''-u{RA) ^ 5].9(r - RA)e'^''' -''>''-u{r) 



Ra 



Y,9{r^RA)e^^^'-^y^- 



Ra 



Ra 

u{r) = 0, 



UBir) ^^g{r- Rb)u{Rb) 



Rb 



E.g(r - RB)u{r) = 



Rb 



Y.g{r-RB) 



Re 



u{r) ~ u{r) = UA{r), 



Ub 



sir) = J2g{r HB)e'(^'-^)-«-^.(HB) ^ Y.sir - RB)e'^''' u{r) 



Rb 



Rb 



Rb 

u{r) = 0. 
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Here we have used first (exploiting tlie hypothesis that the external potential is a very 
smooth function in conparison with g{r)) the property (144) and then (for the quantities 
inside the square brackets) the properties (141) and (143) of the function g{r). In this 
last case the effect of the external potential on the Hamiltonian matrix is only to sum 
the same quantity, u(r), to all the diagonal elements of the matrix, as expected from the 
k ■ p theory (see eq. (106), where the external potential was assumed slowly variable) 



(185) 



u{r) 



7(k^ - iky) 



E 



-f{kx + iky) 


u(r) 














- (r) - 




Ff (r) 




Ffir) 




. Ffir) _ 








u(r) 

7(Ka; - iky) 






-fik.j: + iky) 
uir) 



F^ir) 
F^ir) 
Ffir) 
Ffir) 



Let us note that by reordering the elements of the envelope function vector, we can 
rewrite this equation in the form 



(186) 



E 



uir) 



7(k^ + iky 
Ff (r) • 
Ffir) 
Ffir) 
F^ir) . 





it(r) 
7(k^ - iky) 






7(k^ + iky) 
uir) 




7(K:r - iky 




u{r) 



F^ir) 



F 



A 

K' I 



A (r) 
(r) 

L F^ir) 



which can be more compactly written as 



(187) 



u(r)/ 7(T • k 
7cr • a u(r)/ 



- Ff (r) - 




- Ff (r) - 


Ffir) 
Ffir) 


= E 


Ffir) 
Ffir) 


. F^ir) _ 




. F^ir) _ 



(where / is the 2x2 identity matrix and cr is the vector having as components the 
Pauli spin matrices and ay (44)). This equation is analitically equivalent to the Dirac 
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equation for massless particles (Weyl's equation) of relativistic quantum mechanics (■^); 
therefore eq. (185) is commonly called the Dirac equation for graphene. Since charge 
carriers in graphene obey a relation identical to that describing the relativistic behavior 
of elementary massless spin-(l/2) particles, transport in graphene exhibits many phe- 
nomena, such as Klein's tunneling, analogous to those predicted in relativistic quantum 
mechanics [57-61]. 

Note that in the presence of a magnetic field the operator k = — iV which appears 
in the equation has to be replaced by — iV -I- eA/h, as we have shown in the general 
introduction on the k ■ p method. 

In the absence of an external potential, the quantities ua, u^, ub and u'g arc null 
and thus the matrix equation becomes 



(188) 





-f{kx + iky) 



(r) 

Ff{r) 
Ff{r) 



7(k^ - iky) 










7(/«a; - iky) 






j{kx + iky) 




(r) 
F^{r) 
Ff{r) 
Ffir) 



Since in this case the part of equation corresponding to the point K is decoupled from 
that corresponding to the point K' , we can consider the two parts separately. 
In particular, the part of equation corresponding to the point K is 



(189) 






-f{kx - iky) 




" Ff (r) ■ 


^ E 


' F^{r) ' 


j{kx + iky) 







F§{r) 


F^{r) 



or (using the Pauli spin matrices (44)) 

(190) 7(K,a, + ky<jy)F^{r) ^ -f{k ■ (T)F^{r) = EF^{r). 

This k p Hamiltonian matrix, converted into the momentum representation (see eq. (Ill)), 
has as eigenvalues the dispersion relations of the two degenerate energy bands 
E^ {k) and as eigenvectors the corresponding electron envelope functions F^{r). 
In particular, if we set 



(191) 



det 



7(K:r - iHy) 

7(k^ + iKy) 





' 1 " 




- E 




1- 




1 



(^) For example, compare this equation with eq. (3.62) of ref. [96], with 



0, A 



0, 



eAo = u{r), c substituted by vf = 7/ft, i/'a substituted by [F^ , F^ ]^, and substituted 
by [Fr, Ff]^. 
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we find the dispersion relations 

(192) Ef{K) = s^^kI + «2 = sj\kI 

where s can assume the values +1 or —1. 
If we define the angle a in such a way that 



(193) 
and thus 
(194) 



we have that the corresponding envelope functions (properly normalized, as we will see) 
are 



(195) 
with 
(196) 



Ff^[r) = ^e-'-e'^=('^)i?(-a(K))|s) 



where f2 is the considered surface area, 0s (t) is an arbitrary phase factor and R{a) is a 
spin-rotation operator, given by 



(197) R{a) = 
This can be easily verified noting that 

(198) 7 



e'- 

e-*i 



If^iy 





1 





Kx ~\- il^y 

-i\K\e-''°' 
ikle'" 





1 



/2n 
1 

1 



e"'-'^e"^=i?(-a(M))|s) 
e''^-'-e''^'R{~a{K))\s) 



1 

2\/n 



i\K\e 



2^ 




e~*t 


1 


—is 




e't 


71 


1 


) 
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and also 



(199) E^Ff-^ir)^s^\K\ 
, , 1 



2n 



e-*t 


1 


—is 




e't 


72 


1 


) 



SJ\K\ 



1 



2V^ 



2vn 



76 



— ise 



2V^ 



(where we have used the fact that = (±1)^ = 1 ). 

Instead, the part of equation corresponding to the point K' is 



(200) 



7(^3. + iky) 

^{kx - iky) 



Ff{r) 
Ff{r) 



Ff{r) 
F^'ir) 



or equivalently (using the Pauh spin matrices (44)) 



(201) 7(K,a, - kyay)F^ (r) = 7 





/ 








— Ky 










CT I F^ (r) = EF"^ (r) 



K' 



1 
1 



= 0, 



If we move to the momentum representation (see cq. (Ill)) and enforce 
(202) det|[ , ^ . l{^. + ^^y) 1 _E 

we find the dispersion relations 
(203) 



where s can assume the values +1 or —1. 
The corresponding envelope functions are 



(204) 



1 



/2n 



e"'"-e*('^)i?(a(K))|s), 



with (/>s(/«) an arbitrary phase factor and 



(205) 



IS 

1 



This result is easily verified in a way completely analogous to eqs. (198)-(199) [93]. 

From these functions F^ , Fg , F^ and Fg , we can find the functions ipA and V'b 
and thus the electron wave function ip in the absence of an external potential, using the 
relations (134) and (120). 
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We notice that the energy dispersion relations that we have found in this way near 
K and K' are identical to those one can obtain by first computing the dispersion rela- 
tions in the absence of an external potential by using the nearest-neighbor tight-binding 
technique, and then expanding them near the extrema points. 

Let us now find an expression for the probability density and for the probability 
current density in graphene. 

The probability to find an electron in a region of area S is equal to 



(fi is the area of the whole graphene sheet), where we have exploited the fact that 
each atomic orbital if has a non-zero overlap only with itself (since we use Lowdin 
orthonormalized atomic orbitals) and has significant values only near the atom on which 
it is centered. If the atomic orbital ip is normalized according to (138), the integral 
of its square modulus on $7 is equal to the unit cell area Hq (otherwise, if the usual 
normalization for ip is adopted, this integral is equal to 1 and the following results just 
have to be divided by the constant Hq). Therefore, in this case we have that 





(207) 




Using the relations (134), we have that 



(208) 






Ra Ra 




Fr{RA)Ff{RA) 
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and that 
(209) 



J2\mrb)\^ = Y,rB[RB)^^B{RB) = 

Rb Rb 
Rb 

■\ie'^'e'^-^^F^{RB) + e'^'-^^Ff{RB)\] = 
Y,\F^{Rs)\' + Y.\Ff{Rs)? 



Re 



Re 



''Y,[e^^^'-^y^-Fr{RB)Fr{RB) 

Rb 
Rb 



However the terms contaming the phase factors e^^^ -k) Ra ^ ^i(K ~k)-Rb ^ their com- 
plex conjugates are ncghgiblc with respect to the others. 

Indeed, using the smoothing function g{r), we know from the property (141) with 
r — Ra that X^Jt' 9{Ra — R^a) = 1- Therefore we can insert this sum into the term 



(210) 
obtaining 

(211) 



Ra 



E 



Y,9{Ra-R!a) 



R'. 



e^^^'-^y^-FriRA)Ff{RA) 



that can be rewritten, as a result of the point-symmetry of the function g with respect 
to its center and thus of the fact that g{RA — R'a) = 9{^{Ra — Ra))j this way: 



(212) 



EE5(«A - fi^)e^(^'-^'-«-^Ff *(i2A)Fr (i?^). 

Ra R'. 



If then we use the property (144) with r = R'a and in particular the fact that 

(213) giR'A - RA)FriRA)FfiRA) = giR'A - RA)Fr iR'A)Ff (R'a) 

(due to the smoothness of the envelope functions), the term becomes 



(214) 



Y,[Y,g{RA-RA)e 



i{K'-K)-RA 



^A {R'a)Pa (R'a) 



and, by way of the property (143) with r = R'a, we conclude that the quantities between 
square brackets, and thus the overall term, are very small. 
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Analogously, we can see that the terms 



Ra 



and 

Rb 

are negligible [93]. Since g{r) has non-negligible values only within a few lattice constants 
from its center, the previous considerations are approximately valid also if we limit the 
sums to the atoms contained in the area S. 
We conclude that 



(215) 



f \ij{r)\^dr^no ^ |Va(-Ra)|' + 17o Y1 \^b{Rb 

Ra^S Rb^S 

no E \F^{RA)\' + no \FfiRA)\' 

Raes RAes 
+^0 E \F^iRB)\' + no E \FfiRB)\'o, 

Rb^S Rb^S 

[\F^{r)\^ + (r)p + |Ff (r)p + (r) 



where we have exploited the fact that the envelope functions F are smooth functions, 
which are nearly constant over a unit cell. Therefore we can consider 



(216) 



P = \F^ (r)r + \Ff{rr + \F§ (r)|^ + (r)| 



as a probability density, and the correct normalization condition is 

(217) ^ (li^i'(r)l^ + \Ff{r)\^ + \F^ {r)\^ + \Ff{r)\^) dr = 1. 

We now follow a procedure similar to that used in relativistic quantum mechanics [96] 
to find the expression of the probability current density. Let us consider the envelope 
function equation in the case of long-range external potential (eq. (185)), writing explic- 
itly the operators Kj, and ky (see eq. (155)). Let us consider the time-dependent wave 
function tp{r,t) and thus the time-dependent envelope functions F{r,t) {F will be the 
column vector [F^ , F^, F^ , Fg ]-'"). We now convert the time-independent enve- 
lope function equation into a time-dependent envelope function equation, substituting 
in the r.h.s. of eq.(185) the quantity EF{r) with ih{dF{r,t)/dt) (for stationary states 
^{r,t) = tj)[r)e~^^*'/^ ^ F{r,t) = F{r)e~"^^^/^ , and thus the time-dependent equation is 
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clearly equivalent to the time- independent one). Therefore we can write 



(218) 








OX 


dy 
















4 

ax 


_d_ 

dy 



















^ B 

















-i# + 

OX 


d 
dy 
















ax 


d_ 

dy 









L r g J 



-u{r) 



^A 
^ B 
^A 
^B 



^A 

pK 

^ B 

^A 
^B 



Dividing by 7 and using the Pauli matrices (44) , we can rewrite the equation in this form 
(in the following we will indicate with / the 2x2 identity matrix): 



(219) 



i h 

7 



-i Ux 

—i Ux 

I 
/ 



d_. 
dx 







dt 



u{r) 





/ 
/ 



—1 

dy 



F = 



that, if we define 



iux 

iCTx 



(220) A = 

we can rewrite in this compact way: 



B = 



ioy 



ih f d 



1 \dt 



u{r) 



F = 0. 



7 



If wc left-multiply this equation by the row vector F^ (the conjugate transpose of F) 
we obtain: 



(222) 



F^F = 0. 



Instead, if we consider the conjugate transpose of eq. (219) we obtain 
(223) 



ih f d 



ial 
i crj. 



7 V^* 



/ 
/ 



dy 



^al 
-i at 



I 
/ 



0, 
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which, since crj = ax and crj = (Jy, is equal to 



= 0. 

7 



If we right-multiply this equation by the column vector F, we obtain 



(225) ( l-Ft") AF+( ^F^) BF + —( ^fA F + ^F^F - 

\dx ) \dy J J \dt J 7 

Subtracting (222) from (225), we find 
(226) 



= 0. 



^fM AF + F^'A ( ^F 
ox I \dx 



^FM BF + F^B ( -^F 

dy ) \dy 



+ - 



i h 



7 [\dt 



^ fAf + f^ (If 



dt 







-^(F^AF) + -^(F^BF) + —^{F^F) = 0. 
ax ay 7 at 

Since F^F — P (probability density), we have that (defining vp — l/K) 



(227) -^P = -§i{F^F) = -t (J) V • [{F^AF)x + {F^BF)y 



i VF F^AF)x + {-i VF F^BF)y 



V J, 



which is the well-known continuity equation, if we define as probability current density 
the vector 



(228) 







Jy 





-ivF F'^AF 
-ivF F^BF 



In particular, we have that 
(229) Jx = -ivFF^AF = 



-I Vf 



'I Vf 



^A 



^B 



F 



^B ^A 



^B 



K'* T^K' 



i 

1 
i 
i 



^Ff 



FrFf 



K 



*Ff) 





r -I 

^A 




^ B 




^A 




pK' 
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and that 



(230) Jy 



ivF F^BF = 



-I VF 



-I Vp 



-I Vp 



F 



K* 



^ B 



F^F^ 



F 



K pK' 



F^-F^ 



F. 



K'* 



F, 



F, 



K' 

B 



1 





" 




r 

^ A 


1 










^B 








-1 




^A 





1 







L ^B 






^ B 
-^A 
-^B 

^A 



We note that a different ordering of the elements inside the envelope function vec- 
tor is often used [62,61], in which, instead of F, the vector F = \F^{r), F^{r), 
Fb '{''') ^ ^a" !^)]"^ is considered. Consequently, the k ■ p equation in the case of long- 
range external potential (185) can be rewritten in this way: 



(231) 



w(r) 7(k^ - iky) 
^{kx + iky) u(r) 



E 






F^{r) 
F^{r) 
Ff{r) 
Ffir) 






u{r) 
j{k.j: + iky) 






7(k^ - iky) 
u{r) 



F^{r) 
F§{r) 



B (r) 

A 



Ff{r) 



which is the so-called valley-isotropic representation of the Dirac equation, characterized 
by two identical 2x2 submatrices corresponding to the two valleys K and K' . 

Following this representation, the previously obtained expressions for the probability 
current density can be compactly restated in this form: 



J = VF F [I ® (t)F, 



where I ®(t is the Kronecker product between the 2x2 identity matrix / and the vector 
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(T of Pauli matrices. Indeed, the resulting x and y components of J are 



(232) J^^ vpF {I®(T^)F 



T^K* T^K* T?K'* T^K' 
^B ^A 






1 





" 




r 

^A 


1 













^B 











1 




^ B 








1 







L ^A 



VF ( FrF^ + FrFf + F^Ff + FrF'^' 



B : 



Jy = vpF {I®(ry)F 



T^K* T^K* tpK' 
^A ^B ^B 



^A 






—i 








i 




















—i 








i 








r pK -1 

^A 




^B 




pK' 

^ B 




L ^A J 



-I Vp 



F^-F^F^ + F^ F^ 



which coincide with eqs. (229)-(230). 

It is useful to notice that the Dirac equation in the absence of an external potential 
is not satisfied only by the eigenvector F{r) = [Ff (r), F^{r), F^'{r), F^' [r)]'^ 
with eigenvalue E (as we see in (188)), but is satisfied also by the eigenvector Fi{r) = 
[F^{r), —F^{r), F^ (r), —Fg (r)] ^ with eigenvalue —E, since (188) is equivalent to 








^{k^ - iky) 










- Ff(r) - 


(233) 


-f{kx + iky) 













7(kj, + iky) 




-F^{r) 
Ff{r) 










7(k^ - iky) 







. -Ff{r) _ 



- F^{r) - 

-Ff(r) 
Ffir) 
. -Ff{r) _ 



The wave functions V'('") and ^/'i(r) corresponding to the envelope functions F(r) and 
Fi{r) therefore have opposite energies and thus, being (see eq. (121)) eigenfunctions of 
the Hermitian operator H corresponding to different eigenvalues, are orthogonal. But, 
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due to the form of F{r) and Fi{r) and to eq. (181), we see that ip{r) and have 
the same ipAif) but opposite ipBir). Therefore, if we write the orthogonahty relation 
between 'ij){r) and 7/'i(r), wc have that 

(234) = / ij{r)*i!i{r)dr = 



I \Y.i'A[RAMr -Ra) + Y^ ^BiRsMr - Rb) 

Ra Rb 



dr 



■ [ ^ i;A{RA)^{r -Ra)-Y, i^siRsMr - Rb) 

Ra Rb 

J2\i'A{RA)\' [ W{r - RA)?dr -Y,\^b{Rb)? f W{r- 

R,A R-B ^ 

^ \i^A{RA)?^^ - Hb{Rb)?^0 ^ 



RB)?dr 



Ra Rb 

noY,\i^A{RA)? ^nnY.\'l'B{RB)\^ 
Ra 



Re 



where we have exploited the fact that each atomic wave function ip has a non-zero overlap 
only with itself and has been normalized according to (138). Since (as we have seen) 



(235) 



17o ^ |Va(-Ra)|' ^ f^o 5] |F^(i2A)P + f^o 51 l^f (^^ 

Ra Ra 



Ra 



F^{r)\' + \F^\r)\']dr 



K' 



Qo E \MRb)\' ^ f^o E \FBiRB)\' + E IFfiRs)? 

Rb Rb Rb 

l^{\F^ir)\' + \Ff{rr)dr, 



we conclude that 



(236) / (K(r)p + K'(r)p)dr= / {\F^ {r^ + \Ff {r)\') dr 

and this means that in the absence of an external potential the normalization (217) is 
equivalent to 



(237) 



F^{r)\^ + \Ff{r)Adr 



1 

2' 



^ j^^{\F^{r)? + \F^{r)?)dr^\ 



(the expressions of the envelope functions previously written for graphcnc in the absence 
of an external potential satisfy this normalization criterion). 



546 



P. MARCONCINI and M. MACUCCI 



5. Application of the k ■ p method to carbon nanotubes 



A single-wall carbon nanotube can be described as a graphite sheet rolled, along one 
of its lattice translational vectors (the vector Ch shown in fig. 7), into a cylindrical 
shape [83]. In particular, it is completely specified by the so-called chiral vector Ch, 
which corresponds to a section of the nanotube perpendicular to the nanotube axis and 
thus has a length equal to the nanotube circumference and connects two points of the 
graphene sheet which coincide in the nanotube. This vector can be expressed as a linear 
combination of the real space unit vectors of graphene with integer coefficients n and m 



(238) 



Ch = nai 



ma-) = na 



" V3 - 




' ^/3 - 




2 




2 




1 


+ ma 


1 


= a 


2 




^2 















V3 



[n + m) 



1. 







The corresponding carbon nanotube will be indicated as {n,m). 

If we define the chiral angle of the nanotube (with — 7r/6 < < 7r/6, due to the 
hexagonal symmetry of graphene lattice) as the angle (positive in the clockwise direction) 
between ai and Ch (sec fig. 7) or, equivalently, as the tilt angle of the edges of the 
hexagons constituting the graphene sheet with respect to the direction of the nanotube 
axis, such an angle can be found from the values of n and m noting that 



(239) 



cos 9 



Ch ■ ai 



2n 



|C';i||ai| 2Vn^ + m'^ + nm 



and 



(240) 



sin 9 



{Ch X oi) ■ z 
|C„||ai| 



73^ 



where the right-hand reference frame S' = {x\ y\ z') is that already used in the calcula- 
tions on graphene. In the successive expressions we will identify the previously introduced 
angle 9' with 9' = (tt/G) — 9 (the angle between Ch and the axis x'), as shown in fig. 7, 
and thus we will take the axis x along Ch- 

Following Ando's approach [55,56], the dispersion relations and the electron wave 
functions of a carbon nanotube can be obtained from those of graphene, enforcing for the 
electron wave function the following periodic boundary condition in the circumferential 
direction: 



(241) ^(r + Ch) - i'{r) 

(in the calculations we will not consider the curvature effects (^)). Remembering that 



(*) For the effects of tlie finite curvature on the electronic properties of carbon nanotubes see, 
for example, ref. [97] and the references therein. 
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using the tight-binding technique the electron wave function can be expressed as 
(242) V(r-) = i'AiRAMr ~Ra) + Y. i^siRsMr - Rb), 



Ra 



Rb 



the boundary condition can be written as 



(243) 



V'(r + CO = 

^ ^A{RAM{r + Ch) -Ra) + J2 ^B{RB)^{{r + Ch) - Rb) = 



Re 



R 



i^AiRAMr - {Ra - Ch)) + Y i'B{RB)^{r - {Rb - C,)) 
J2i^A{{RA - Ch) + Ch)^{r ~ {Ra ~ O) 
+Y,i^B{{RB - Ch) + Ch)^{r - {Rb - Ch))= 

Rb 

Y^a{R*a + Ch)v{r -R*a) + J2 '^b{R*b + Ch)^{r - R*b) = 

R\ R% 

^{r) ^VM(i?:i)^(r - R*a) + YMR*B)v>{r - R*s) 



R% 



(where we have used the fact that, being C/i a hnear combination with integer coefficients 
of the real space lattice unit vectors, also Ra — Ch and Rb ^ Ch ai'c atomic positions, 
defined R*a and -R^)- Thus the boundary condition is equivalent to the two conditions 



(244) 



^A{RA + Ch) = ^A{R*A), 
^B{R*B + Ch) = ^B{R*B)- 



If we use the expressions (134) for iPa{i') and 'ipB{f) (and we define again the generic 
atomic position Ra and Rb, instead of R\ and R*b), these conditions can be rewritten 
in the following form: 



(245) 



f e'^-(«^+^'')Ff (il^ + Ch) - i e^'^'e'^'-(^-^+^'')Ff' (i2^ + C„) = 
e'^-^^F^{RA) -ie'^'e'^'-^^Ff{RA), 
ie'^'e'"^ ^^F^{Rb + Ch) + e»^'-(«-+^")Ff' (J?b + Ch) = 

• ^^K■RB pK ( ) _^ g.K'RB pK' ^j^^y 

Multiplying the first equation of (245) by g{r — RA)e^''^''^^ , summing it over Ra and 
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then using the properties of the function^ (defined in eqs. (f4f ), (f 43) and (f 44)), we find 

(246) e'^^-^-^5(r-i2A)Ff (il^ + CO 
Ra 

Ra 

^ gir ~ Ra)F^{Ra) - ^ e^"' ^ 9{r ~ i?A)e'(^'-^)-«-<' {Ra) =^ 
Ra 



Ra 



JKCh 



Ha 



n{r + Cu) 



Ra 



5]5(r--i2A)e'(^'-^)-«- 



Ra 



Ff{r + Ch) 



^g(r-fi^)e*(^'-^)-^- 
Ra 



F^\r) 



e'''-''-F^{r + Ch)=F^{r). 
If we calculate the scalar product between K and Cu we obtain 

(247) 



3 ^ ' 3 ' 



where m ~ n ~ ZN + i^, with = or ±1 and N a proper integer. Therefore we have 
that 



(248) 



and thus the first boundary condition near K is 

(249) e»"i"Ff (r + C,) = Ff(r), 
or equivalently 

(250) F^{r + Ch)^e-''^F^{r). 

Multiplying the second equation of (245) by g{r — RB){~ie~'^^ e~^^'^"), summing it 
over Rb and then using the properties of the function t/, we find analogously [93] 

(251) e^^-^'-Ff (r + C„)=^f (r). 

Substituting the value of e^^'^'^ ^ we can rewrite this boundary condition in the form 

(252) e'"5"^^f (r + C,) = Ff (r), 



THE k ■ p METHOD AND ITS APPLICATION TO GRAPHENE-RELATED MATERIALS 



549 



or, equivalently 

(253) F^{r + C,) = e-''^F^{r) 

Thus the periodic boundary condition near K is 



(254) 



F^ir + Cu) 





" F^{r) - 







which can be written in this compact way: 
(255) F^(r + C,0 = 



F^(r). 



However, as we have previously seen (eq. (195)), in the absence of an external potential 
the envelope functions have the following form: 



(256) 



1 



'2U 



R{-a{K))\s) 



/2M 



with the surface area VI = LI, where L = \Ch\ and £ is the length of the nanotube. Thus 
the periodic boundary condition becomes 

(257) ^^e*'^-(''+<=''')e*"^-'^'^)i?(-a(K))|s) = e''^ ^=e"^-'' e"^'^"^^ R{~a{K))\s) , 
V2M V2l£ 

or equivalently 
(258) 

This condition can be written also in the following way 
(259) 



3 e 



z27rn 



or, equivalently 
(260) 
and thus 
(261) 



with fi integer. 

This condition on Kx can be obtained also in a different way, enforcing the boundary 
condition on the overall wave vector k. In order to do this, we have to observe that, 
considering only the periodic lattice potential inside the graphene sheet, the wave function 



27rzy 

K^L = h 2iin 

o 
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ipir) has to be a Bloch function r)e*'" '', where u{k,r) has the periodicity of the 
lattice. 

Thus the boundary condition 

(262) V(r + Ch) = Hr) 
is equivalent to 

(263) u{k, r + Ch)e'''-'-''+'=''^^ = u{k, r)e*'="'. 

Since we know that u{k,r) has the lattice periodicity and thus u{k,r + Ch) = u{k,r) 
{Ch being a linear combination with integer coefficients of the lattice unit vectors) the 
boundary condition can also be written as 

(264) 6**=-'^'' = 1, 
or, equivalently 

(265) kCh = 27rTO. 

Thus the boundary condition is (being Ch = Ch/\Ch\ = Ch/L) 

and (using eq. (247)) 

(oai^ 27r, 2^_ K Ch 2^. 2^ ~ 2^ 
zo7 K,„ = — m — (K r = — m = — m iV v = 

27r / _ ~ u\ 27r / _ v\ , , 

(with fi = m — N) , which is equal to the previously found expression. 

If we substitute this condition on k^. in the dispersion relations of graphene, we find 



(268) Efnii^v) = si\i^\ s-i^ kI + kI = s-i^ Ky{hY + k: 



where s = +1 and s = — 1 indicate the conduction and valence bands, respectively. 

We notice that now ky is the wave vector k of the nanotube, which, being a substan- 
tially unidiniensional material, has a one-dimensional Brillouin zone with width 27r/T 
(where T is the length of the unit cell of the nanotube, along its axis, which can be easily 
found from the numbers n and m characterizing the nanotube [83]). Correspondingly, 
Ky is the difference between the wave vector k of the nanotube and the component of K 
along y. 
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As to the envelope functions near K, if, starting from eq. (195), we choose as value 
of the arbitrary phase (f>s = —a/2 and then we enforce the condition on Kx, we can write 



(269) Ff^ir) 



1 



/2Le 





1 






/n 




1 


■ gg-»(f+a) " 


2^ 


/n 


1 




1 


sb^{n, Ky) 


2i 


/n 


1 



e * 2 





1 


—is 











1 












a) 




2\ 









-ise 
1 



The function bv{n,Ky) = e *(2+") can be found noting that a has been defined (see 
eq. (193)) in such a way that 



(270) 

this means that 
(271) 



Hx + iKy = |K|e'(2"''"); 



i(f +q) _ + lHy 



and thus 

(272) 6,(n, Ky) = e-<^+") = (e'(^+")) " 



We can proceed analogously for the boundary conditions near if'. 
Indeed, multiplying the first equation of (245) by g{i — JRyi)(ie^*^ e~^^ '^'^), summing 
it over Ra and then using the properties of the function g, we find [93] 

(273) e*^'-^''<'(r + C„ ) = Pf (r). 

The scalar product between K' and Ch is equal to 



(274) 



where we have used the previously introduced relation m — n ~ 3N + ly with = or 
±1 and TV a proper integer. Thus we have that 
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and consequently the boundary condition near K' is 

(275) e-*^F5(r + C,0 = F^'(r), 
or, equivalently 

(276) F^'{r + Ch) = e^'"^ Pf {r). 

On the other hand, multiplying the second equation of (245) by g{r — Rb)^'^^ ^ 
summing it over Rb and then using the properties of the function g, wc find [93] 



(277) 



Substituting the value of e^^ ''^'^ we can rewrite this second boundary condition near 
K in the form 



(278) e-'^Ff{r + Ch)^F§'{r) 
or, equivalently 

(279) Ff' (r + CO = e'^^^f (r). 
Thus the overall periodic boundary condition near K' is 



K' I 



(280) 



F^'{r + Cn) 
which can be written in a compact form 



(281) 



[r + Ch) = e'T"F^ (r) 



- ttK' , 



Substituting the form that, in the absence of an external potential, the envelope functions 
have near K' (eq. (204)) 



(282) Ff(r) 



ml 



e^'^-'^e*('^)i?(a(K))|s) 



ml 



i?(a(K))|5), 



the periodic boundary condition becomes 
1 



(283) 



V2M 

or, equivalently 
(284) 

This can be rewritten in the form 



R{a{K))\s) 



mi 
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or, equivalently 

27rj/ 

(286) K^L = — + 27rn 

O 

and thus 

(287) K^ = ^(n+'£j ^ k^{n), 
with n integer. 

Analogously to what we have done near K, this condition on Kx can be found [93] 
also setting e*'" '^'' = 1. 

If we substitute this condition on Kj. in the dispersion relations of graphene, we find 

(288) Ef^nii^y) = sj\k\ = s^^kI + kI = s-i^k,(nY + nl, 

where ky now is the wave vector k of the nanotube and Ky is the difference between the 
wave vector k of the nanotube and the component of K along y. 

On the other hand, if, starting from eq. (204), we choose as arbitrary phase (j)s = a/2 
and then we enforce the condition on , we find [93] as envelope functions in the carbon 
nanotube near K' 



(289) Ff^ (r) = 



Sby(n, Ky) 

1 



where (using the definition of the angle a: see eq. (193)) 



(290) 6,(n, Ky) = e<f +") = "" f + = 



ky (n) 



\J + K^ \Jky{n)'^ + K^ 



If jTi — 71 is a multiple of 3 and thus = 0, for ri = and rl = we have that Ky{h) = 
and ki,{n) ~ 0, and consequently — s7|kj,| , which vanishes for Ky ~ 0, so that 
= E- = 0. This means that when to — rt is a multiple of 3 the points K and K\ 
where the upper and lower bands of graphene are degenerate, are among the values of k 
allowed by the periodic boundary condition, and thus the nanotube is metallic. 

Instead, if to — n is not a multiple of 3 and thus v ~ ±1, the allowed fcs nearest 
to K and K' correspond to n = and n = 0, for which Kv{h) = ^2tt/{3L) and 
ky{n) ~ ±2tt/{3L), and consequently 



(291) E,^sj\l{^y + Kl. 

In particular, the minimum and maximum values of the nanotube bands are obtained 
with the further position Ky = and therefore are equal to 

(292) E, = sj^; 
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Fig. 8. - The nanotube (10,0) and its dispersion relations, obtained both by means of the tight- 
binding method (sohd hnes) and (for the bands corresponding to the smallest values of |K^(n)| 
and |Kiy(n)|) by means of the k ■ p method (dashed lines). 

10 

8 

■k 6 

w , 
O 4 
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2 



-3 -2-10123 

E/Yo 

Fig. 9. - The density of states per unit length of the nanotube (10,0), obtained both by means 
of the tight-binding method (solid lines) and (in a smaller region around i5 = 0) by means of 
the k ■ p method (dashed lines). 




thus the bandgap of the nanotube is 



(293) E, = E^-E^ = 2j'=^ = ^ = = 2^^jo "^"^ 

^ ' 9 + '3L 3L 3L 2 L y/3 



where dt = L/tt is the nanotube diameter. Therefore we have that the bandgap of the 
nanotube depends on the reciprocal nanotube diameter. 

We can observe that the approximate approach for the computation of the density 
of states in carbon nanotubes proposed by J. W. Mintmire and C. T. White [98], being 
based on a linear approximation of the dispersion relations of graphene near the extrcma 
points, can be seen as a consequence of a fc • p study of the nanotube energy bands. 

In fig. 8 we compare the dispersion relations that we have obtained for the same 
carbon nanotube using the nearest-neighbor tight-binding method and the k ■ p method 
(without considering curvature effects) [27,99-101]. We see that the k p method gives a 
good approximation for the portions of energy bands of the nanotube deriving from the 
graphene dispersion relations around K and K' . 

In fig. 9, instead, for the same nanotube we show both the density of states that we 
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have obtained by properly differentiating the tight-binding dispersion relations, and the 
density of states deriving from the Mintmire- White approach [27,99]. We see that this 
last approximation gives good results near E = 0, thus in the region corresponding to 
the graphene dispersion relations around K and K' . 

6. Application of the k ■ p method to graphene nanoribbons 

A graphene sheet can be laterally confined (along the y-direction) to form a graphene 
nanoribbon (extending in the x-direction) . The properties of the nanorribbon strongly 
depend on the characteristics of the boundary. Here we will consider nanoribbons with 
perfect zigzag and armchair edges, that can be easily studied using the Dirac equation and 
enforcing the correct boundary conditions [63,64,102-106]. An analysis of the boundary 
conditions that have to be enforced in nanoribbons with more general terminations can 
be found in ref. [107]. In particular, we will perform the analytical calculations in the 
absence of an external potential following Brey and Fertig's approach [63,64], but using 
the representation adopted in the previous sections. While inside the nanoribbon each 
atom has 3 nearest-neighbor atoms, for the atoms on the edges of the ribbon some of 
the nearest-neighbor lattice sites are outside the ribbon and thus are not occupied by 
a carbon atom. These lattice sites are instead occupied by passivation atoms (such as 
hydrogen atoms), which saturate the dangling bonds. The correct boundary condition 
to be enforced in our calculations is the vanishing of the wave function in correspondence 
of these lattice sites (let us call them "boundary lattice sites" ) . 

6'1. Zigzag nanoribbons. - In the case of zigzag nanoribbons (fig. 10), the graphene 
sheet has been cut at an angle of 30° with respect to the nearest-neighbor carbon bonds, 
and therefore the edges have a zigzag shape. In order to simplify the following calcula- 
tions, we can choose (sec fig. 10) the graphene lattice vectors in the real space ai and 02 
(and consequently those in the reciprocal space bi and 62) in this way (we express them 



A A A A A A A 




a 



Fig. 10. - Sfcetch of a zigzag nanoribbon with zigzag lines (the black atoms are carbon atoms, 
while the grey atoms are passivation atoms). 
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in the reference frame S = {x, y, £)): 



(294) 



ai 



bi 



2 

V3 

a 

2 



2n 
a 

2-K 




02 



ho = 



2 

V3 

a 

2 



27r 
a 

2tt 




(which, being hi = 27r(a2 x i)/(ai • {a^ x £)) and 62 = 27r(£ x ai)/(ai • (02 x £)), fulfih 
the relation • bj — iudij). Consequently we have that 



(295) 



1 47r 
K = -{b2-bi) = - 
6 s 3a 



1 Att 
K' = Uhi-h,)^^ 
6 s 3a 



" -1 " 




" -K ' 


















" 1 " 




" K ' 












1 













where we have defined K =^ 47r/(3a). For our choice of ai and a2, the angle 9' from the 
vector fli + 02 {i.e. from the axis x used in previous calculations) to the axis x (taken 
in the longitudinal direction) is equal to tt/2. 

Therefore the total wave function is given by (eq. (120)) 
(296) ^{r) = ^A{RA)v{r - Ra) + Y. ^siRsMr - Rb). 

Ra Rb 
with (eq. (181) with 9' = 7r/2) 



V;A(r) = e'^-'^i^f (r)+e^^-"-Fj^-(r), 



iK' r t:^K' I 



(297) 



s'-^'-Ff (r) + e'-^''-i^f' (r), 



where (using eq. (295)), if we write r = [x,y,0] we have that K ■ r = —Kx and that 
K' ■ r ~ Kx. In the absence of an external potential, the envelope functions satisfy the 
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usual Dirac equation (eq. (188)) 



(298) 










OX 


_9_ 

dy 
















7 


ax 


a 

' dy 



















F^{r) 

















~i # 4 

ax 


_d_ 

dy 




(r) 














ox 


_d_ 

dy 












\ F^[r) 1 
Ff(r) 
(r) 
- Ff{r) _ 



Due to the translational invariance along the a;-direction, we can write the envelope 
functions as the product of a propagating part along the longitudinal direction x and of 
a confined part along the transverse direction y. Therefore we can assume that 



(299) 



Air) 

K 

B 



and that 



F^ (r) 
Ff (r) 



n (y) 

•ff (y) 



We have to enforce that the overall wave function vanishes in correspondence with the 
"boundary lattice sites" on the lower and upper edges of the ribbon. Let us define as 
W the real width of the nanoribbon, i. e. the distance between the lowest row of carbon 
atoms (all of type A) and the highest row of carbon atoms (all of type B); if the ribbon 
has N zigzag lines across its width, we have that W = {3N — 2)ac-c/2- If we take 
y = in correspondence of the row of "boundary lattice sites" on the lower edge, the 
row of "boundary lattice sites" on the upper edge will be for y = Vt^ = + 2ac-c = 
(3A^ + 2)ac_c/2. The proper boundary condition thus implies that, for every x, ?/'(x, y = 
0) = i/'(x, y = W) ~ 0. Since in the zigzag nanoribbon all the "boundary lattice sites" on 
the lower edge belong to the B sublattice, while all those on the upper edge belong to the 
A sublattice, looking at eq. (296) and observing that the atomic orbitals cp are strongly 
localized around the atom on which they are centered, the boundary condition on the 
wave function is equivalent to setting, for every x, TpB^x^y = 0) = il)Aix,y = W) ~ 0. 
Using eq. (297), we have that 

(300) 7/'B(a;,?; = 0) = Vx -e-*^^F|^(2;, ?/ = 0) + e'^^Fj^' (x, y = 0) = 

_^~rKx^rK,x^K ^q) ^ ^rKx ^ik'^x ^K' (q) = Q Vx ^ 

<I>f(0) = 0, $f'(0)-0 



and that 



(301) i;A{x,y ^W) Vx^ e-'^^'^ F^ {x,y = W) + e'^'' Ff {x,y = W) = 
e-'^'=e"'^''^^{W) + e*-^^e*''-^$f' (ly) = Va; ^ 
$f(VF) = 0, $f'(VF)=0. 
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As we can see, in zigzag nanoribbons the boundary conditions do not couple the envelope 
functions relative to the Dirac points K and K' . 

First let us make the calculation around the point K. The corresponding part of the 
Dirac equation is: 



(302) 



dx 



dx 



d_ 

dy 



d_ 

dy 



dx 



dx 



^ n r 

dy 



^ n r 

dy 



(r) 
Ff (r) 



= E 



Ff(r) 
Ff(r) 



(y)e^«^- + e*'===-^$f (y) 



dy 



-f 

dy 



E 



- (r) " 


= E 











d 

dy 



4t. 



E 

1 



^f{y) 



which can be rewritten as 



(303) 



dyj 



7 
7 



Obtaining $§^(y) from the second of (303) and then substituting ^^{y) from the first 
of (303), we find: 



(304, = 2 („. + A) *f fa, . (2)' („. + A) - A) *f (,) = 
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the solution of which is 



(305) $g(y) = yle^y + Be-^«, 



^X2 



with z=yK^— ( — j ( and thus E ~ :^'^\/ k\ — 



Substituting (y) back into the first of (303), we obtain that 



(306) $f(y)^2(^^^_^^<i>f(y) = 

l{{K,-z)Ae'-y + {ti,+z)Be-'-y). 



Let us now enforce the boundary conditions on ^g{y) and (y) 



(307) $f (0) = 0^ A + B = O^B = -A; 

>^^{W) = 0^1 ((k, - z)Ae/'^ + (k, + z)Be-''^) = ^ 

(k, - z)Ae'^ - (k, + z)Ae-'^ = ^ 
(k^ - z)Ae''^ = {k^ + z)Ae-''^ ^ 

g-2zM/ _ l^x — Z 
~t~ Z 



As we can see, in zigzag nanoribbons the longitudinal and the transverse wave vectors 
are coupled. 

Incidentally, note that, instead of eq. (307), an equivalent equation can be used [106]; 
indeed, being E ~ ±7\/k^ — z'^ and thus {E/"/)"^ = i^x ^ -^^j have that 



(308) 



e 
E 



2zW 



z {k^-z){k^ + z) 



±{K, + z)e-'^. 



{kx + zY + zf 



Here we consider real values of k^. 

If we graphically represent (fig. 11) the two functions f\{z) = e^^^^ and ji{z^ = 
(kx — z)I(kx + 2), we see that (apart from z = 0, which corresponds to identically null 
<I>'s) there is an intersection between /i and ]i for a real value of z (and thus eq. (307) 
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-1 











Fig. 11. - Graphical solution (in the real domain) of eq. (307) (the dotted lines are the asymptotes 

0f/2(^)). 



has a real solution z) only if k^; > and if fi{z) is steeper than /2(z) hi z = 0, i.e. if 



(309) 



dz 



h{z) 



> 



z=0 
J z=0 

Kr + Z + Z 



dz 



> 



1 









- 2 = 







Kx + Z {Kx + Zf' 
[Kx + Zf 



- 1 1 
2W > ^ > — ^ Kx> ^. 

Kx W 

If instead < l/W, eq. (307) does not have real solutions z (apart from z = 0). 
In the case of real z, from eq. (307) we can find that 



(310) 



-2zW 



z^K,(l-e-2-^)=z(l + e-'^^^) 
z 



1 + e 



-2zW 



^zW ^ p-zW 



1 _ 6-2214/ g2iy _ ^-zw tanh(zl^) 
(z = does not have to be considered) and thus 

~,2 I ^^rA^2 



(311) 



.2 2 



tanh^(zVF) 
/ cosh^(zW^) - sinh2(zll^)\ 



_^2^^2 / cOsh^W \ ^ 



sinli^(zVl^) 



sinh2(zM^) 



sinh^(zW^) 







z 


7 




sinh(ziy) 



Since (for the properties of the hyperbolic sine function) | sinh(ziy)| > |ziy| = |z|iy, we 
see that in this case 



(312) 



E 

7 



1 



< 



\A 



\z\W W 
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We can write (exploiting what we have found from the boundary conditions) that 

E 

1 
E 

1 



(313) $f (y) - ((«, - z)Ae"^ + («, + z)Be-'y) 
I {{k^ - z)Ae'y - («, + z)Ae~'y) 



[n^ie^y - e-'^y) - z(e^y + e^^^')) = 

— 2A sinh(zj/) ~ z co&h{zy)) = 

1 z \ 

2A — — sinh(zw) — zcosh(zw) I = 

E Uanh(zVF) ^ ^ ^7 

2A ^ cosh{zW) sinh(zy) — sinh(zM^) cosh(zy) 

E^ sinh(zl^) 

2A cos\v{zW) sinh(— zy) + sinh(zM^) cosh(— zy) 

sinh(ziy) 

-2^ {l?—T^lT7^^ sinh(z(I4^ - y)) = 



E sinh(zVK) 

~2Asign('-— f-^) sinh(z(W^ - y)), 
V 7 sinh(zW^) / 

where in the last step we have taken advantage of the fact that, due to eq. (311), the 
product between 7/i? and z/sinh(zW^) can only be equal to +1 (if the two quantities 
have the same sign) or — 1 (if they have opposite signs) . 
Moreover we have that 

(314) $f (y) = Ae'y + Be""^ = Ae'^ - Ae-'-y = A(e"^ - e""^) = 2Asinh(zy). 

These are edge states, each one exponentially localized on one edge of the ribbon. 

These edge states correspond to bands flattened towards E = 0, as we can see 
both from the graphical solution of eq. (307) (where we observe that we have an in- 
tersection between /i and /2 for a z coordinate very close to and thus the energy 
E = ±7-\/k^ — has a very small value), and from our previous analytical conclusion 
that \E/^\ < 1/W in this case. Since the Dirac point K, folded into the Brillouin zone 
(— 7r/a, 7r/a) of the zigzag nanoribbon (the unit cell of which is of length a), corresponds 
to kx = — 47r/(3a) + 27r/a = 27r/(3a), the condition Kx > l/W (under which we have a 
real solution and thus the edge states) is equivalent to kx = Kx + Kx > 27r/(3a) + 1/W 
(note the difference between the total wave vectors k and the wave vectors k measured 
from the Dirac points). Therefore in the region 27r/(3a) + 1/W < kx < Tr/a we have 
two bands flattened towards E ~ 0; this means that the zigzag nanoribbons arc al- 
ways metallic [108]. However, further studies [109-111] have shown that actual zigzag 
nanoribbons have a non-zero gap deriving from a staggered sublattice potential due to 
edge magnetization. 

Let us now instead consider the imaginary solutions (with K„ real) of 

eq. (307). In this case the dispersion relation E = ±7-^/ — becomes E = 
±7-v/ -|- K^, from which we see more clearly that Kx and have the meaning 
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of longitudinal and transverse components of the wave vector, measured from the Dirac 
point. The solutions arc given by 



(315) 



-2zW '^x 



g — i2Z(Ka:+iK„) _ g — i2 Z (k x +iKn ) gi2-;rm 

KnW = 1{Kx + iHn) — T^fn => tan(K„W^) = — 



tan(K„PF) 



(with m integer); k„ = corresponds to identically null $'s and thus docs not have to 
be considered. We have that 



(316) 



2 I 2 



tan( 

cos^(k„H^) + s'm^ (KnW) 



^ V , ^2 _( C0S^(^«W") \ 2 

[^nW)) +"""lsin^(.„W^)+'] 



sin [KnW) " sin2(K„#) 



sin(K„W^) 



since (for the properties of the sin function) | sin(K„M^)| < |k„VK| = |k„|PF, we see that 
in this case 



(317) 



> 



1 



\nn\W W 

We can write (exploiting what we have found from the boundary conditions) that 
(318) $5^(y) = ^ {{n^ - iK„)Ae>^-y + [k^ + zK„)Se-*'^"^) = 



E 

1 
E 

E 

2iA^ I — sin(K„-v) — k„ cos(k„u) 

E Vtan(K„W^) ^ ^ ^' 

7 cos(K„Ty) sin(K„?/) - sin(K„Ty) cos(K„y) 

ZlA — K,, 



((k, - tKn)Ae'''"y - {k^ + iKn)Ae-'^"y) 



2iA {k^ sin(K„?/) - k„ cos(K„y)) 



-2iA 



7 f^r, 



~2iA sign 



E sm{KnW) 

E Kr, 



7 sin(K„H/) 



sin(K„(M^ - y)) = 
sin(K„(W^-y)), 



where in the last step we have taken advantage of the fact that, due to eq. (316), the 
product between j/E and K„/sin(K„M^) can only be equal to +1 (if the two quantities 
have the same sign) or — 1 (if they have opposite signs) . 
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Moreover we have that 



(319) 



$f (y) = Ae'^-y + Be- 



= At 



These are clearly confined states extending all over the ribbon. 

The calculations around the point K' are completely analogous. The corresponding 
part of the Dirac equation is 



(320) 



dx 



dx 



d_ 

dy 



d_ 

dy 



dx 



dx 



d_ 

dy 



d_ 

dy 







(r) 
Ff (r) 



= E 



Ff{r) 
Ff{r) 



4$f (y)e*'=-^ + e* 



E 



l^x dy 



/ d_ 

X dy 



Ak^x d 



d 

dy 



dy 



d^n (y) 
" $f (y) 



E 



(r) 



E 

7 



it' 



(y) 



which can be rewritten as 



(321) 



-.+^'1 $r(y) = ^$f(y), 



7 

4-7^)$?' (y) = -<(y)- 



Obtaining (y) from the second of (321) and then substituting (y) from the first 
of (321), we find 



(322) c&f (y) 



E 



dy 



d 



d'- 



dy 



E 



dy dy^ 



(y)=(-) $f (y), 



F 



dy 

<ff' (y) 



$r(y) 
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the solution of which is 

(323) $f' (y) = Ce/-'y + Der'-'y , 



with z' = \l k'^'^ ^ ^ ^ ^^'^ ^^^'^ ^ ^ ^'x^ 



Substituting [y) back into the first of (321), we obtain that 



(324) Ci,^ (y) = ^|^4 + ^^$f'(2;)^ 

^ ( n'Ce^'y + n'Oe-'-'y + z'Ce'-'y - z'De-''y 
E \ 



1_ 

E 



((4 + z')Ce^'y + (4 - z')^e-^'^) . 

Let us now enforce the boundary conditions on (y) and (y): 
(325) $f' (0) = ^ C + Z? = ^ D = -C: 



$f (M/) = ^ -1 ((4 + z')Ce^ + (4 - ^')^e-^ ^ j = : 

(4 + z')Ce"''^ - (4 - z')Ce-"'^ = ^ 

(4 + z')Ce''^ = (4 - z')Cer''^ ^ 
4 + z' (-4) - z' 



(-4) + z" 



which is equal to eq. (307) if we substitute with — k^. Therefore the calculations are 
completely analogous to those seen around the point K . 
We consider again real values of k'^. 

We conclude [93] that (apart from z' = 0, which corresponds to identically null $'s) 
eq. (325) has a real solution z' only if —4 > i-e- if < ~i/W. 

If instead k'^ > ~1/W, eq. (325) docs not have real solutions z' (apart from z' = 0). 
In the case of real z', from eq. (325) we can find that [93] 



(326) 



tanh(zW) 



{z' = does not have to be considered) and thus [93] 
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The corresponding $ functions are [93] 

(328) $f (y) = + z')Ce^''' + (4 ^ z')De-^'y) 



2C sign 



E 



7 sinh(zW) 



smh{z' {W - y)); 



$f (y) =Ce''y + De-' y ^ 2Csinh(z'y). 

These are edge states, each one exponentially localized on one edge of the ribbon. 

Also in this case, these edge states correspond to bands flattened towards E = 0. Since 
the Dirac point K' , folded into the Brillouin zone (— 7r/a,7r/a) of the zigzag nanoribbon, 
corresponds to — Att / {3a) — 2tt / a = — 27r/(3a), the condition k'^ < —1/W is equivalent 
to k'^ = + k'^ < — 27r/(3a) — 1/W. Therefore also in the region —ir/a < k^ < 
— 27r/(3a) — 1 /W we have two bands flattened towards E ^ 0, which confirms the metallic 
nature of zigzag nanoribbons. 

Let us now instead consider the imaginary solutions z' — in'^ (with k'^ real) of 

eq. (325). The dispersion relation E = ±7-^ k'J^ — 2'^ becomes E = ±7^ 
The solutions are given by [93] 



(329) 



tan«M^) 



(kJj = corresponds to identically null $s and thus does not have to be considered) and 
thus [93] 



(330) 



/ 2 J 2 



sin^«ty) 



sin«M/) 



> 



1 



Wjw w 



The corresponding $ functions arc [93] 

(331) c&f iv) = I ((4 + *<)Ce-> + (4 - i<)De— 



2iC sign 



E k' 



7 sin{K'^W) 



sin«(M^-y)); 



(y) Ce"'-y + De-'^'^y = C2i sinin'^y) . 

These are confined states extending all over the ribbon. 

Obviously, once the expressions of the functions $ have been obtained, the overall 
wave function is given by the equations (296), (297) and (299). 

In fig. 12 we show the bands of a zigzag nanoribbon with = 45 zigzag lines and 
of a zigzag nanoribbon with = 50 zigzag lines, that we have computed both with a 
simple tight-binding model not including edge magnetization effects (thick dotted lines) 
and with the k ■ p (Dirac equation) method (thin solid lines). For low energy values and 
for not too narrow ribbons the results obtained with the two techniques are very similar. 
In both cases, the presence of the two bands flattened towards zero and corresponding 
to the edge states can be clearly seen. 
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k,(A-') k.(A-') 

Fig. 12. - Bands of a zigzag nanoribbon with A'^ = 45 zizag lines (a) and with A'^ = 50 zigzag 
lines (b), computed both with a simple tight-binding model not including edge magnetization 
effects (thick dotted lines) and with the k ■ p method (thin solid lines). The two dashed lines 
correspond to the energy values iLj/W; the dispersion relations in the region between the two 
dashed lines are obtained for real values of z, while those outside this region correspond to 
purely imaginary values of z. 



6'2. Armchair nanoribbons. - Instead, in the case of armchair nanoribbons (fig. 13), 
the graphene sheet has been cut along the direction of the nearest-neighbor carbon bonds, 
and therefore the edges have an armchair shape. In order to simplify the following 
calculations, we can choose (see fig. 13) the graphene lattice vectors in the real space 
ai and a2 (and consequently those in the reciprocal space bi and 62) in this way (we 
express them in the reference frame S = {x,y,z)): 
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Fig. 13. - Sketch of an armchair nanoribbon with dimer lines (the black atoms are carbon 
atoms, while the grey atoms are passivation atoms). 
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(which, being bi = 27r(a2 x z)/{ai ■ (02 x £)) and 62 = 27r(£ x ai)/(ai • (02 x £)), fulfih 
the relation • bj = 2T:5ij). Consequently we have that 























(333) 
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4(5. 


-61) 


_ 47r 
5] 3a 


-1 




-K 
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K' 




-b2) 


47r 
5] 3a 


1 





K 



1 



with K = 47r/(3a). For our choice of ai and 02, the angle 9' from the vector ai + 02 (i.e. 
from the axis x used in previous calculations) to the axis x (taken in the longitudinal 
direction) is equal to 0. 

Therefore the total wave function is given by (eq. 120) 



(334) 



V'(r) 



E 

Ra 



Rb 



with (eq. (181) with 6' = 0) 

r ^A{r) = e'-f^-Ff (r) - i e^^'^Ff' (r), 



(335) 



V;B(r) = ze^-^-Ff (r) + e^^'-F|^'(r), 



where (using eq. (333)), if we write r = [a;,?/,0]"'" wc have that K ■ r = —Ky and that 

K' ■ r = Ky. In the absence of an external potential the envelope functions satisfy the 
usual Dirac equation (eq. (188)) 



(336) 7 







dx 



_§_ 

dy 



E 






F^{r) 
F^{r) 
(r) 



d_ 

dy 







_§_ _ 

dx 



d_ 

dy 






_d_ _ 

dx 



_§_ 

dy 



Ff(r) 
Ff (r) 
Fj^'(r) 



Due to the translational invariance along the x-direction, we can write the envelope 
functions as the product of a propagating part along the longitudinal direction x and 
of a confined part along the transverse direction y. Here we have to consider the same 
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longitudinal component for the wave vector measured from K and K' because in this 
case if we consider k'^ ^ Kx the boundary conditions are satisfied for every x only by the 
identically null wave function. Therefore we can assume that 



(337) 



and that 



(r) 
Ff (r) 



We have to enforce that the overall wave function vanishes in correspondence with the 
"boundary lattice sites" on the lower and upper edges of the ribbon. Let us define as 
W the real width of the nanoribbon, i.e. the distance between the bottom row and the 
top row of carbon atoms of the ribbon; if the ribbon has N dimer lines across the ribbon 
width, we have that W = {N — l)a/2. If we take y = in correspondence of the row of 
"boundary lattice sites" on the lower edge of the ribbon, the row of "boundary lattice 
sites" on the upper edge of the ribbon will he aX y = W — W 2a/2 = W + a = 
(N + l)a/2. Therefore, for every x, we must have ^{x^y = 0) = ^'(x,y = W) = 0. 
We notice that in an armchair nanoribbon the "boundary lattice sites" on the lower and 
upper edges belong to both the A and the B sublattices. Therefore, looking at eq. (334) 
and observing that the atomic orbitals if are strongly localized around the atom on which 
they are centered, the boundary condition on the wave function is equivalent to setting, 
for every x, tpA{x,y = 0) = -tpBix.y = 0) = ipAix.y = W) = ipBix,y = W) = 0. Using 
eq. (335) we obtain the following 4 boundary conditions: 



(338) iI;a{x, y = 0) = Vx ^ e^'^^Ff (x, y = 0) - ie'^°Ff{x, y = 0) = 

"'-''^f{0) = Vx 



(x, y = 0) - zFf (x,y = 0) = e^-==-$f (0) 



(0) 



zc&f (0) = 0; 



(339) i^Bix, y = 0) = ^x ^ ie 



iKO T7K 



B 



(x,y = 0) + e^^OFf' (x,y 



zFf (x, y = 0) + (x, y = 0) = ze---<I>g (0) 



z$f(0) + <i>f (0) = 0; 



0) = 

$5 (0) = Vx: 



(340) i^A{x,y = W) = Q Vx 



-iKW ipK 



(x, y = W)- ze'^^Ff' (x, y = W) 



-iKW iK^x 



(W) = Vx 



^^{W)-ie'^^<pf{W) = 



(341) iPBix,y = W) ^0 Vx=^ie-*^'^F^(x,y 



-iKW ^iK^x jf.K 



iKW ^iKxX jf.K 



le 



-■W) + 
Vx 



e^^^Ff'(x,y = #) 



As wc can see, in armchair nanoribbons the boundary conditions couple the envelope 
functions relative to the Dirac points K and K' . 
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We can solve the part of the Dirac equation around the point K, that is 



(342) 



dx 



_d_ 

dy 



_a d_ 

dx dy 





F^{r) 
F^{r) 



E 



Ff (r) 



repeating the calculations made for zigzag nanoribbons (eqs. (302)- (306)) and obtaining 
that 



(343) 



^Aiy) = I - ^)Ae^'' + i^x + z)Be-^y) 
(y) = Ae/y + Be-'^y , 



with z = \l k1 — ^-^^ and thus E = ±7-\/k2 — z^. 

Analogously, we can solve the part of the Dirac equation around the point K' , that 



(344) 



ox 



dx 



d_ 

dy 



d_ 

dy 



Ff{r) 
Ff{r) 



= E 



Ff{r) 
Ff{r) 



repeating the calculations made for zigzag nanoribbons (eqs. (320)-(324), with the dif- 
ference that and z' here have to be replaced by Kx and z) and obtaining that 



(345) 



(y) - ^ {{^^ + ^)Ce''' + {^^ - z)De-^y) , 

*f' (y) = Ce'-y + Der'-y , 



with (as written before) z ~ \j Kj.^ — vT j thus E = k^^ — z^. 

Let us define z = iKn- In this case the dispersion relation becomes E ~ 
±7V + Kn^; therefore Kx and k„ = —iz are the longitudinal and transverse com- 
ponents of the wave vector, measured from the Dirac point. 

The functions $ become 



(346) 



' $f (y) = {{i^x - in^)Ae'^-y + {kx + iK„)Be-"^''^) , 



$f (y) = Ae'^'-^y + Be^'^-y 



I {y) ^ Ce'^-y + De 
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Now we can enforce the 4 boundary conditions (338)-(341), obtaining: 

(347) $A (0) - * *A ' (0) = ^ 

{Kx - iKn)A + {Kx + iK,n)B - i{Hx + iHn)C - i{Kx - iKn)D = ^ 

Kx{A + B-iC - iD) + iKn {-A + B - iC + iD) = 0; 

(348) (0) + $f' (0) ^Q^i{A + B) + {C + D) = Q^ 
A + B-iC -iD = Q- 

(349) e-'^^<^f{W) - ie'^'^^fiW) = =^ 

e-'^^{nx - in^)Ae"'-'^ + e''^^ {n^ + in,,)Be-'^-^ 

(350) ie-'^^$f (ty) + e'^"'$f' (ly) = 0^ 

■g-iifiy^^g»«„M/ Be-'^"^) + e*^^(Ce*«"^ + De"*''"^) = ^ 

In the following we examine the different cases in which all of these 4 boundary conditions 
are satisfied. 
Case I 

If K„ = the condition (347) is equivalent to the condition (348), and the condition (349) 
is equivalent to the condition (350). 
But the condition (348) is satisfied if 



(351) A + B -iC -iD = 0^ A + B = iC + iD ^ 



A + B = G 
C + D = -iG 



(where wc have defined A + B = G). 

The condition (350) instead is satisfied if (exploiting the fact that k„ = 0) 

(352) AeK'^"-^)^ + iJe-»(-"+^)i^ - zCe^(-"+^)^ - zZ?e-^('="-^)^ - ^ 

{A + B)e~''^^' - i{G + D)e:^^ = ^ Ge"*^^ - Ge''^^ = ^ 
-G(e'^^' - e-^^^) = ^ -G2i sin(W) = 0. 



THE k ■ p METHOD AND ITS APPLICATION TO GRAPHENE-RELATED MATERIALS 



571 



Since in this case (k„ = 0) for G = all the $ functions (346) would become identically 
null and thus we have to consider G 0, this equation can be satisfied only if sm{KW) = 
0. But, since K 47r/(3a) and W ^ {N + l)a/2, we have that 



(353) sin(i^W) = ^ sin ( ^{N + f )| j = ^ sin ( y (iV + 1) ) = 



and this is true only if TV + 1 is a multiple of 3, i.e. if iV + 1 = 3AI (with M integer) and 
thus N = 3M — 1. In this case we have that 



(354) 



2tt , , 27r, 

KW = —{N + l) = —i'SM) = 27rM =^ 

K = 2mZ- 2mZ- ~ K = 0(= Kn) 
WW ^ ' 



and the nanoribbon is metallic (as we will see). Being k„ = 0, the $ functions (346) are 
equal to 



(355) 



^Aiv) = ^{^xA + K^B) = ^K^{A + B) = ^K,G, 
$f(2y)=^ + B = G, 

^a\v) = ^(«:xC + K.,D) = ^Ac,(C + D) = -^K,tG, 
[ <i,fiy)^C + D^-iG. 



Case II 

The other possibility is to satisfy the conditions (347)-(348) in this way: 



(356) 



A + B-iC -iD = j 2B-2iC = j = 

-A + B -iC + iD = {) \ 2A-2iD = {) \ D = -iA 



(where in the first step we have summed and subtracted the two equations of the system), 
and to satisfy the conditions (349)- (350) enforcing 



(357) 



Using (356), we can write these equations in the following form: 
(358) <^ 
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If now we separate the real and imaginary part of k„ (i.e. we write k„ as k„ = K„,. + iK„ 
we have that 



(359) < 



_QQ-K.„iW ^i(k„^+K)W _|_ j^^K„iW^-i(K„,-K)W ^ 

[^e-K„,w cos((K„r - K)W) + Be'"-^ cos((k„^ + K)W) 
_5e-K„.H' cos((K„r + K)W) - Ae«"'"' cos((k„^ - K)W)\ 
+i [Aer^^^^ sin((K„^ - K)W) - Be''"'^ sin((K„^ + K)W) 



-Be 



sin((K„r + K)W) + Ae 



sm 



((k„, - K)W)] = 0, 



[ - ^e"'"-^ cos((k„^ - K)W) + Se"^-^ cos((K„r + K)W) 
_5e-«-^ cos((k„^ + + Ae'^'"^ cos((k„^ - is:)!^)] 

+i [ - Ae-""'"^ sm{{Knr - K)W) ~ Se"-"' sin((K„r + K)W) 
_5e-K„.vi/ sin((K„^ + K)W) - Ae«"''^ sin((K„^ - K)W)] = 0, 



^'"^) [Acos((K„r - K)W) ~~Bcos{{K„r + K)W)] 



+ e" 



)[Asin((K,„. - - Bsin((K„^ + K)W)] = 0, 

(e''-^ - e-"-'^) [A cos((K„r -~K)W) + B cos((K„r + 
_j(e'=n.H' g-K„.w') sin((K„^ - if + B sin((K„^ + K)W)] = 0, 

-2 sinh(K„jM^) [A cos((K„r - K)W) ~ B cos((K„r + K)W^ 
+i2cosh(K™W^)[Asin((K„r - if)!!") - B sin((K„r + iC)W')] = 0, 
2sinh(K„,W^) [Acos((ft„r - K)W) ^ B CQ^[{Kr,r ^ K)W)\ 
^ -i2cosh(K„jM^)[Asin((K„^ - K)W) + Bsin((K„^ + K)W)\ = 0. 

If wc sum and subtract the two equations, we obtain 



(360) < 



' 4:Smh{KmW)B cos{{K„r + K)W) 

i4cosh{K„iW)B sin{{Knr + K)W) = 0, 
4sinh(K„iI^)Acos((K„r - K)W) 
+i4cosh(K„iW^)Asin((K„r - K)W) = 0, 

_B[sinh(K„,W^) cos((K„r + K)W) — i cosh(K„,W^) sin((K„r + = 0, 

A[sinh(K„iM^) cos((/v„r. — K)W) — i cosh(K„iM^) sin((K„r — A')W^)] = 0. 



Apart from the case A = B = 0, which (being also C = —iB and D = —iA) gives 
identically null functions $, both of these two equations are satisfied in 3 cases, that we 
will indicate with II-A, II-B and II-C. 
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(361) 



siuh{KmW) cos((K„r + K)W) - i cosh(K„jT^) sin((K,ir + K)W) = 0, 



[ Smh{KniW) cos((K„r — K)W) — i COSh{KniW) sin((K„r — K)W) = 0. 
If we separately equate to zero tlie real and imaginary parts, we find 



(362) 



smh{KmW) cos((K„r + K)W) — 0, 
cosh{KniW) sin((K„r + K)W) = 0, 
sinh{KmW) cos((K„r — K)W) — 0, 
cosh^KniW) sin((K„r — K)W) = 0. 



Since the hyperbolic cosine is never equal to zero, these become 



(363) 



' sinh(K„,W^) cos((k„^ + K)W) = 0, 
sm{(Knr + K)W) = 0, 
sinh(K„iT^) cos((K„r — K)W) = 0, 
L sin((K„, - K)W) = 0. 



However, when the sine of an angle is equal to zero, the cosine of that angle is certainly 
different from zero; therefore the previous equations become 



(364) 



sinh(K„iW^) = 0, 
sin((K„^ + K)W) = 0, 
sm{{K„r - K)W) = 0. 



Since the hyperbolic sine is null only when its argument is null, we conclude that in this 
case: 



(365) 



am{{K„r + K)W) = 0, 
sin((K„, - K)W) = 0, 



sm{{K„ + K)W) = 0, 
sin((K„ - K)W) = 0. 



From the condition on sin((K„ + K)W) it follows that 

(366) sin((K„ + K)W) = ^ (k„ + K)W = titt 
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(with n integer). Then from the condition on sin((K,i — K)W)^ substituting what we 
have just found and then remembering that K = 47r/(3a) and that W = [N + l)a/2, we 
obtain that 



(367) sin((K„ - K)W) = ^ sin ( ( n-^ - K - K ]W ] = fim[m: - 2KW 



, 4:71 , , a 

sm mr - 2 — {N + 1)- 
3a 2 



n — 4 IS an integer. 



sin TT n — 4 



1 



This is true only if iV + 1 is a muhiple of 3, i.e. if + 1 = 3M (with AI integer), i.e. 
ii N = 3M — 1; this means that the nanoribbon is metalhc (as we will see). In this case 
the $ functions (346) are equal to 



(368) 



^Aiy) = I {{^. - in^)Ae'^-y + (k, + inr,)Be-'^-y) , 
$f (y) = Ae*'*"^ +56-'"-^, 



*f (2/) = -h ii^- + ^>^n)Ce'^-y + [k, ~ iKn)De-'^-y) 



$f' (y) = Ce'^-y + De-'^-y = -i (Be"'"y + Ae-"""'^ 



7 



E 

-i^ {{kx + iKn)Be''^'-y + (k^ - iK^)Ae~'''^"y) 



that can be written as a superposition of the modes 



(369) { 





^(^x + iK„)yle 








i«:„)Se-'""^ 


1>f(2/) = 








= Be-'^'^y 




-^(k^ +m„)iAe*'^"^ 


and < 


7 . 
- ^(A^. 










+ iKn)iBe'''-y 


.*f' (j/) = 


-iAe'^'-^y, 




(2/) 







with K„ = {mr/W) — if and k„ = — k„. We notice that, since in this case N + 1 = 3M, 
we have that 



(370) 



KW 



4tt 



-{N + l) 



a 2tt 



3a' '2 3 
and therefore k„ can be written as 



2-K 

(A^ + 1) = — (3M) = 2-kM => K = 2M^ 



W 



(371) 



nZ- - K 
W 



TT 



w 



2K - K 



W 



W 



-n-^ + 4Af-^ ] - K= (4M - n)-^ -K = n-^ - K, 



W 



with fi = AAI — n integer. Clearly, if k„ satisfies E — ±7\/K^JJ^~r«Vi^, also k„ = — k„ 
satisfies = ±7-\/k^^^"T"k^. 
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It can be observed that in the particular case in which k„ = we find again Case I. 
Case II-B 

Equations (360) are satisfied also if 



(372) 



smh{KmW) cos((K„r + K)W) — i cosh{KmW) sin((K„r + K)W) = 0, 



If we separately equate to zero the real and imaginary parts of the first equation, we find 



(373) 



s'mh{KniW) cos((K„r + K)W) = 0, 
cosh{KniW) sin((K„r + K)W) — 0, 
^ = 0. 



Since the hyperbolic cosine is never equal to zero, these become 



(374) 



sinh(K„,iy) cos((k,„. + K)W) = 0, 
sin((K„r + K)W) = 0, 
^ = 0. 



But when the sine of an angle is equal to zero, surely the cosine of that angle is different 
from zero; therefore the previous equations become 



(375) 



smh{KniW) = 0, 
sm{{Knr + K)W) = 0, 
A = 0, 



Since the hyperbolic sine is null only when its argument is null, we conclude that in this 
case: 



(376) 



Sm{{Knr + K)W) = 0, 

A = 0, 



K„ real, 

sm{{Kn+K)W) = 0, 
A^O. 



Due to the fact that A = 0, also D = -iA = (while C = ~iB). 
Instead the consequence of the condition on sin((K;„ + K)W) is 



(377) 



sin((K„ + K)W) = ^ (k„ + K)W = mi 

K„ + ii = K„ = n— K. 

W W 
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In this case the $ functions (346) are equal to 



(378) 



' (y) = ((«.^ - i«„)Ae"^"^ + («. + z«„)Se-"'"^) 

*f (y) = I ii^x + i^n)Ce'''-y + (k^ - iKn)De-'^"y) 
III 



[ $f' (y) = Ce'^'-y + De-'""" 



-iBe' 



Case II- C 

Finally, eqs. (360) are satisfied also if 



(379) 



B = 0, 



I sinh(K„iM^) cos((K„r — K)W) — icosh(K„iW) sin((K„r — K)W) = 0. 
With calculations analogous to Case II-B, we conclude [93] that in this case: 



(380) 



B = 0, 

sin((K„ - K)W) = 0. 



Due to the fact that B = 0, also C = -iB = (while D = -iA). 
Instead the consequence of the condition on sin((K„ — K)W) is 



(381) 



sin((K„ - K)W) = ^ (k„ - K)W = mr 
K„ — K = K„ = n— — h A. 

w w 



In this case the $ functions (346) are equal to 



(382) 



*A (y) = ^ {{^x - iKn)Ae'''-y + (k. + iKr,)Be-''^-y) = 



(y) = I {{^'sc + in,,)Ce'^"y + {k, - tK,,)De-'^-y) = 

— (kx iKn)iAe "y = p;('*2; -\- ik,n)iAe ~ "y , 
h/ III 



I $f (y) = Ce^'^"^ + Z?e- 



-iAe' 



~iAe'^ 
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with 
(383) 



nZ- + K 
W 



W 



K = h— K 

W 



(where n = —n is an integer). Clearly, if k„ satisfies E = zt'j\/Kx^ 
satisfies E = ±.^^Kx^ + k^- 

In conclusion, in all the cases we have that 



I ^ , also 



(384) 



(?/) = ^(kx +«K„)Ae- 



with A being a proper normalization constant, k„ = [n-K jW^—K and E = ±7%/ k^.^ + k„2. 
Consequently, for eq. (335) we have that 



(385) 



^Air) = e^^-Ff (r) - z e'^'^i^/ (r) = 
(e-'^^<i>^(y) - ze'^'^$f (y)) e'«^-- = 



-i{Kn-{-K)y\ AkxX 



-^{kx + iKn)A2is\n{{Kn + K)y) e*"-'' 



and that 
(386) 



iK r Tj^K 



F#(r) + e'^'-Fj^'(r 



(^ie-'^'^2^$f (y) + e'-f^y$f' (y)) e*'*=="^ = 



-iA I 



-17421 sin ((k„ + K)y) 6***== 
2Asin((K„ +i^)y) 6*'^==^. 



We observe that in large ribbons the lowest-energy modes will have values of k„ much 
smaller than K and thus their wave functions will be characterized by a transverse 
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wave vector approximately equal to K and by a transverse wave length about equal to 
2'k/K = 27r (3a/(47r)) = 3a/2, i.e. of the order of the lattice constant. 
No edge state exists in armchair nanoribbons. 

Using the relations K = Att/ (3a) and W = {N + l)a/2, we have that 

_ TT _ n27r 47r _ 27r(3n - 2(7V + 1)) 

- " " ~ (iV + l)a ~ 3^ ~ 3(7V + l)a " 

Since En = ±7A/Ka;^ + Kn^, we have a double band degeneracy if, for any integer n, 
another integer n' exists such that k„' = — k„ and thus E^' = E^. This happens if 

(388) 3n' - 2{N + 1) = -(3n - 2{N + 1)) 3n' = -3n + 4(7V + 1) ^ 

A{N+1) 



with n and rt' integer, which means that + 1 has to be a multiple of 3, i.e. iV + 1 = 3M 
with M integer, or equivalently N = 3M — 1 (so that n' = —n + AM). 

We also observe that among the allowed k„'s (given by eq. (387)) we have k„ = if 
an integer n exists, such that 

(389) 3n - 2(iV + 1) = ^ n = ^^^^ , 

which again means that + 1 has to be a multiple of 3, i.e. iV + 1 = 3Af with AI integer, 
or equivalently N = 3 A/ — 1 (so that n — 2M). 

Therefore an armchair nanoribbon has a double band degeneracy and has k„ = 
among the allowed values of k„ only if it has a number of dimer lines TV = 3M — 1 (with 
M an integer). In this case for k„ = we have E = ±j\kx\ which vanishes for — > 
and thus the nanoribbon is metallic. Instead for ^ 3AI — 1 the armchair nanoribbon 
is not metallic and has non-degenerate bands. 

This conclusion is coherent with the fact that the dispersion relations of an armchair 
nanoribbon can be obtained from those of graphene enforcing the Dirichlet boundary 
conditions at y = and y = W; this means that there has to be an integer number of 
transverse half-wavelengths Ay/2 inside W; thus it must happen that 

A„ 27r TT 

(390) w^n^^ky = — =n-^ 

Z Ay W 

(where ky is the transverse component of the total wave vector, measured from the origin 
of the reciprocal space). Therefore the bands of the ribbon can be obtained by cross- 
sectioning those of graphene along the parallel lines ky = mr/W, and then folding them 
into the Brillouin zone (— 7r/(A/3a), 7r/(-\/3a)) of the armchair nanoribbon (the unit cell 
of which has a length 3ac-c ~ V^a). There are bands of the nanoribbon with a zero 
gap, and thus the nanoribbon is metallic, only if some of the lines with ky = nn/W pass 
through a Dirac point of graphene (where the graphene dispersion relations have a zero 
gap). But, since 

rrr , '^W AtT 47r iV + 1 „ TV + 1 TT 

(391) W^iN + l)-^a=——^K= — = ——^=2- 



2 N + 1 3a 3 2VF i W 
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Fig. 14. - Bands of an armchair nanoribbon witii = 98 dimer lines (a) and with A'' = 99 dimer 
Unes (b), computed both with a tight- binding method not including the reduction of the bond 
lengths at the edges (thick dotted lines) and with the k ■ p method (thin solid lines) . 



this is possible only if -I- 1 is a multiple of 3, i.e. -|- 1 = 3M with Al integer, or 
equivalcntly N = 3M — 1. 

A more exact tight-binding analysis (taking into consideration the reduction of the 
carbon-carbon bond lengths parallel to dimer lines at the edges with respect to the bond 
lengths in the core of the ribbon) leads to the appearance of a small gap also in this 
subset of armchair nanoribbons, that have to be more correctly considered as quasi- 
metallic ones [109-111]. 

In fig. 14 we show the bands of an armchair nanoribbon with A^ = 98 dimer lines 
(metallic) and of an armchair nanoribbon with A^ = 99 dimer lines (semiconductor), 
that we have computed both with a tight-binding method not including the reduction 
of the bond lengths at the edges (thick dotted lines) and with the k ■ p (Dirac equation) 
method (thin solid lines). As we see, for low energy values and for not too narrow ribbons 
they are nearly coincident. 

All previous considerations are valid both for real values of (propagating modes), 
and for purely imaginary values of (evanescent modes). 

As an application of the relations (229) and (230) to the case of an armchair nanorib- 
bon in the absence of an external potential, we can observe, using the (384) (with k„ 
real and real or purely imaginary), that 

(392) J. = V, [prF^ + FrF^ + FrFf + vr ) = 

+{k* - iK„)«A*e-*'""^e~"^>(-i)Ae"'"«e*'^-^ 

Vp ^|Ape*^''""'*-^^(K* - IKn + + IKn + K* - «K„ + + iK^) = 
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which if Kx is real (and thus k* = k^) is equal to (remembering that vp = ^/h) 

2 

(393) J.='ivF^\A\^K,^A\A\^^K,, 

while if Kx is purely imaginary (and thus k* = —k^) is null. 

Note that (using (216) and (384)) if Kx is real the probability density is equal to 



2 



(394) P = \F^{r)\^ + \F^ (r)p + |F|^(r)p + |Ff (r 

2 / 'v \ 2 



(I) |«;.+*A^„|2|A|2 + |A|2+(2) \Kx+tK,,\'\A\^ + \Af' 

2\A\'- (l+(lY\^x+ ^^nA = 2\A\' ( 1 + ^'^^it^"^ 



^2 ~ 



2|A|Ml + |!) =2|Ap2 = 4|Ap 



Moreover, since in this case the energy dispersion relations are E = ±j\/kx'-^ + k„2, the 
mean velocity of the electrons is 

1 dE 1 / 1 1 

Therefore if k^. is real we have that = Pvx, as expected. 

As to the transversal part of the probability current density, we have that 



(396) J, = ^ VP [FrF^ - ^^f - Ff ' Ff' + F^ ' Ff' 
-i t;F ^ (« - iK„)A*e''''"''e-<^^e-*'""2^e*'"="^ 

- iK„)a*e-*'""?'e-*«-^(-i)Ae'«"2'e"=-^ 
+a*e~"^"^e-<^(K:, +iK„)(-i)Ae*''"^e"'-^) = 

-i VF ^ e''^''""''-'^(K* - - Ka; - iK„ " K* + ZK„ + K^; + ZK„) = 0, 

III 

as expected in a transversally confined structure. 



7. Conclusion 



The fe • p method and the related envelope function method are widely used to study 
the physical properties of materials within a continuum approach, without having to re- 
sort to an atomistic analysis, which requires (in the case of large structures) a prohibitive 
computational effort. These methods have been developed in many and sometimes quite 
different ways by several authors and have been successfully applied to a multitude of 
different problems. This explains the great variety and inhomogeneity of the related 
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literature. In this review, we have briefly described the basics of these methodologies, 
dwelling upon the treatments that we have considered more useful for an easy compre- 
hension. For a detailed explanation of the different approaches, the interested reader can 
resort to the many papers and books on the topic, some of which are listed in the refer- 
ences. In particular, we have focused on the application of the k ■ p method to graphcnc 
and graphene-related materials, where it results in a description of the electronic proper- 
ties in terms of the Dirac equation. We have compared the different formulations adopted 
in the literature and we have shown how this continuum approach allows to quickly ob- 
tain the most relevant electrical properties of graphene, carbon nanotubes and graphene 
nanoribbons. 
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